text:
@testset "Decision Model kwargs" begin
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys5_re = PSB.build_system(PSITestSystems, "c_sys5_re")

    @test_throws MethodError DecisionModel(template, c_sys5; bad_kwarg=10)

    model = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

    model = DecisionModel(
        MockOperationProblem,
        get_thermal_dispatch_template_network(
            NetworkModel(CopperPlatePowerModel; use_slacks=true),
        ),
        c_sys5_re;
        optimizer=GLPK_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    model = DecisionModel(
        get_thermal_dispatch_template_network(),
        c_sys5;
        optimizer=GLPK_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

    #"Test passing custom JuMP model"
    my_model = JuMP.Model()
    my_model.ext[:PSI_Testing] = 1
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(
        get_thermal_dispatch_template_network(),
        c_sys5,
        my_model;
        optimizer=GLPK_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    @test haskey(PSI.get_optimization_container(model).JuMPmodel.ext, :PSI_Testing)
    @test (:ParameterJuMP in keys(PSI.get_optimization_container(model).JuMPmodel.ext)) ==
          false
end

@testset "Set optimizer at solve call" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    template = get_thermal_standard_uc_template()
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "test"),
    )
    UC = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    output_dir = mktempdir(cleanup=true)
    @test build!(UC; output_dir=output_dir) == PSI.BuildStatus.BUILT
    @test solve!(UC; optimizer=GLPK_optimizer) == RunStatus.SUCCESSFUL
    res = ProblemResults(UC)
    @test isapprox(get_objective_value(res), 340000.0; atol=100000.0)
    vars = res.variable_values
    @test PSI.VariableKey(ActivePowerVariable, PSY.ThermalStandard) in keys(vars)
    @test size(read_variable(res, "StartVariable__ThermalStandard")) == (24, 6)
    @test size(read_parameter(res, "ActivePowerTimeSeriesParameter__PowerLoad")) == (24, 4)
    @test size(read_expression(res, "ProductionCostExpression__ThermalStandard")) == (24, 6)
    @test size(read_aux_variable(res, "TimeDurationOn__ThermalStandard")) == (24, 6)
    @test length(read_variables(res)) == 4
    @test length(read_parameters(res)) == 1
    @test length(read_duals(res)) == 0
    @test length(read_expressions(res)) == 1
    @test read_variables(res, ["StartVariable__ThermalStandard"])["StartVariable__ThermalStandard"] ==
          read_variable(res, "StartVariable__ThermalStandard")
    @test read_variables(res, [(StartVariable, ThermalStandard)])["StartVariable__ThermalStandard"] ==
          read_variable(res, StartVariable, ThermalStandard)
    @test read_parameters(res, ["ActivePowerTimeSeriesParameter__PowerLoad"])["ActivePowerTimeSeriesParameter__PowerLoad"] ==
          read_parameter(res, "ActivePowerTimeSeriesParameter__PowerLoad")
    @test read_parameters(res, [(ActivePowerTimeSeriesParameter, PowerLoad)])["ActivePowerTimeSeriesParameter__PowerLoad"] ==
          read_parameter(res, ActivePowerTimeSeriesParameter, PowerLoad)
    @test read_aux_variables(res, ["TimeDurationOff__ThermalStandard"])["TimeDurationOff__ThermalStandard"] ==
          read_aux_variable(res, "TimeDurationOff__ThermalStandard")
    @test read_aux_variables(res, [(TimeDurationOff, ThermalStandard)])["TimeDurationOff__ThermalStandard"] ==
          read_aux_variable(res, TimeDurationOff, ThermalStandard)
    @test read_expressions(res, ["ProductionCostExpression__ThermalStandard"])["ProductionCostExpression__ThermalStandard"] ==
          read_expression(res, "ProductionCostExpression__ThermalStandard")
    @test read_expressions(res, [(PSI.ProductionCostExpression, ThermalStandard)])["ProductionCostExpression__ThermalStandard"] ==
          read_expression(res, PSI.ProductionCostExpression, ThermalStandard)
    @test length(read_aux_variables(res)) == 2
    @test first(keys(read_aux_variables(res, [(PSI.TimeDurationOff, ThermalStandard)]))) ==
          "TimeDurationOff__ThermalStandard"
    export_results(res)
    results_dir = joinpath(output_dir, "results")
    @test isfile(joinpath(results_dir, "optimizer_stats.csv"))
    variables_dir = joinpath(results_dir, "variables")
    @test isfile(joinpath(variables_dir, "ActivePowerVariable__ThermalStandard.csv"))
end

@testset "Test optimization debugging functions" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    template = get_thermal_standard_uc_template()
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "test"),
    )
    model = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    container = PSI.get_optimization_container(model)
    MOIU.attach_optimizer(container.JuMPmodel)
    constraint_indices = get_all_constraint_index(model)
    for (key, index, moi_index) in constraint_indices
        val1 = get_constraint_index(model, moi_index)
        val2 = container.constraints[key].data[index]
        @test val1 == val2
    end
    @test get_constraint_index(model, length(constraint_indices) + 1) === nothing

    var_keys = PSI.get_all_variable_keys(model)
    var_index = get_all_variable_index(model)
    for (ix, (key, index, moi_index)) in enumerate(var_keys)
        index_tuple = var_index[ix]
        @test index_tuple[1] == PSI.encode_key(key)
        @test index_tuple[2] == index
        @test index_tuple[3] == moi_index
        val1 = get_variable_index(model, moi_index)
        val2 = container.variables[key].data[index]
        @test val1 == val2
    end
    @test get_variable_index(model, length(var_index) + 1) === nothing
end

@testset "Decision Model Solve with Slacks" begin
    c_sys5_re = PSB.build_system(PSITestSystems, "c_sys5_re")
    networks = [StandardPTDFModel, DCPPowerModel, ACPPowerModel]
    for network in networks
        template = get_thermal_dispatch_template_network(
            NetworkModel(network; use_slacks=true, PTDF=PTDF(c_sys5_re)),
        )
        model = DecisionModel(template, c_sys5_re; optimizer=ipopt_optimizer)
        @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        @test solve!(model) == RunStatus.SUCCESSFUL
    end
end

@testset "Default Decisions Constructors" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model_ed =
        EconomicDispatchProblem(c_sys5; output_dir=mktempdir(), optimizer=HiGHS_optimizer)
    moi_tests(model_ed, false, 120, 0, 120, 120, 24, false)
    model_uc =
        UnitCommitmentProblem(c_sys5; output_dir=mktempdir(), optimizer=HiGHS_optimizer)
    moi_tests(model_uc, false, 480, 0, 240, 120, 144, true)
    ED_output =
        run_economic_dispatch(c_sys5; output_dir=mktempdir(), optimizer=HiGHS_optimizer)
    UC_output =
        run_unit_commitment(c_sys5; output_dir=mktempdir(), optimizer=HiGHS_optimizer)
    @test ED_output == RunStatus.SUCCESSFUL
    @test UC_output == RunStatus.SUCCESSFUL
end

@testset "Test Locational Marginal Prices between DC lossless with PowerModels vs StandardPTDFModel" begin
    networks = [DCPPowerModel, StandardPTDFModel]
    sys = PSB.build_system(PSITestSystems, "c_sys5")
    ptdf = PTDF(sys)
    # These are the duals of interest for the test
    dual_constraint = [[NodalBalanceActiveConstraint], [CopperPlateBalanceConstraint]]
    LMPs = []
    for (ix, network) in enumerate(networks)
        template = get_template_dispatch_with_network(
            NetworkModel(network; PTDF=ptdf, duals=dual_constraint[ix]),
        )
        if network == StandardPTDFModel
            set_device_model!(
                template,
                DeviceModel(PSY.Line, PSI.StaticBranch; duals=[NetworkFlowConstraint]),
            )
        end
        model = DecisionModel(template, sys; optimizer=HiGHS_optimizer)
        @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        @test solve!(model) == RunStatus.SUCCESSFUL
        res = ProblemResults(model)

        # These tests require results to be working
        if network == StandardPTDFModel
            push!(LMPs, abs.(psi_ptdf_lmps(res, ptdf)))
        else
            duals = read_dual(res, NodalBalanceActiveConstraint, Bus)
            duals = abs.(duals[:, propertynames(duals) .!== :DateTime])
            push!(LMPs, duals[!, sort(propertynames(duals))])
        end
    end
    @test isapprox(LMPs[1], LMPs[2], atol=100.0)
end

@testset "Test ProblemResults interfaces" begin
    sys = PSB.build_system(PSITestSystems, "c_sys5_re")
    template = get_template_dispatch_with_network(
        NetworkModel(CopperPlatePowerModel; duals=[CopperPlateBalanceConstraint]),
    )
    model = DecisionModel(template, sys; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    @test solve!(model) == RunStatus.SUCCESSFUL

    res = ProblemResults(model)
    container = PSI.get_optimization_container(model)
    constraint_key = PSI.ConstraintKey(CopperPlateBalanceConstraint, PSY.System)
    constraints = PSI.get_constraints(container)[constraint_key]
    dual_results = PSI.read_duals(container)[constraint_key]
    dual_results_read = read_dual(res, constraint_key)
    realized_dual_results =
        read_duals(res, [constraint_key])[PSI.encode_key_as_string(constraint_key)]
    realized_dual_results_string =
        read_duals(res, [PSI.encode_key_as_string(constraint_key)])[PSI.encode_key_as_string(
            constraint_key,
        )]
    @test dual_results ==
          dual_results_read[:, propertynames(dual_results_read) .!= :DateTime] ==
          realized_dual_results[:, propertynames(realized_dual_results) .!= :DateTime] ==
          realized_dual_results_string[
              :,
              propertynames(realized_dual_results_string) .!= :DateTime,
          ]
    for i in axes(constraints)[1]
        dual = JuMP.dual(constraints[i])
        @test isapprox(dual, dual_results[i, :CopperPlateBalanceConstraint__System])
    end

    system = PSI.get_system(model)
    parameter_key = PSI.ParameterKey(ActivePowerTimeSeriesParameter, PSY.PowerLoad)
    param_vals = PSI.read_parameters(container)[parameter_key]
    for load in get_components(PowerLoad, system)
        name = get_name(load)
        vals = get_time_series_values(Deterministic, load, "max_active_power")
        vals = vals .* get_max_active_power(load) * -1.0
        @test all(vals .== param_vals[!, name])
    end

    res = ProblemResults(model)
    @test length(list_variable_names(res)) == 1
    @test length(list_dual_names(res)) == 1
    @test get_model_base_power(res) == 100.0
    @test isa(get_objective_value(res), Float64)
    @test isa(res.variable_values, Dict{PSI.VariableKey, DataFrames.DataFrame})
    @test isa(read_variables(res), Dict{String, DataFrames.DataFrame})
    @test isa(PSI.get_total_cost(res), Float64)
    @test isa(get_optimizer_stats(res), DataFrames.DataFrame)
    @test isa(res.dual_values, Dict{PSI.ConstraintKey, DataFrames.DataFrame})
    @test isa(read_duals(res), Dict{String, DataFrames.DataFrame})
    @test isa(res.parameter_values, Dict{PSI.ParameterKey, DataFrames.DataFrame})
    @test isa(read_parameters(res), Dict{String, DataFrames.DataFrame})
    @test isa(PSI.get_resolution(res), Dates.TimePeriod)
    @test isa(get_system(res), PSY.System)
    @test length(get_timestamps(res)) == 24
end

@testset "Solve DecisionModelModel with auto-build" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    template = get_thermal_standard_uc_template()
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "test"),
    )
    UC = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    output_dir = mktempdir(cleanup=true)
    @test_throws ErrorException solve!(UC)
    @test solve!(UC; optimizer=GLPK_optimizer, output_dir=output_dir) ==
          RunStatus.SUCCESSFUL
end

@testset "Test Serialization, deserialization and write optimizer problem" begin
    fpath = mktempdir(cleanup=true)
    sys = PSB.build_system(PSITestSystems, "c_sys5_re")
    template = get_template_dispatch_with_network(
        NetworkModel(CopperPlatePowerModel; duals=[CopperPlateBalanceConstraint]),
    )
    model = DecisionModel(template, sys; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=fpath) == PSI.BuildStatus.BUILT
    @test solve!(model) == RunStatus.SUCCESSFUL

    file_list = sort!(collect(readdir(fpath)))
    model_name = PSI.get_name(model)
    @test PSI._JUMP_MODEL_FILENAME in file_list
    @test PSI._SERIALIZED_MODEL_FILENAME in file_list
    ED2 = DecisionModel(fpath, HiGHS_optimizer)
    @test build!(ED2, output_dir=fpath) == PSI.BuildStatus.BUILT
    psi_checksolve_test(ED2, [MOI.OPTIMAL], 240000.0, 10000)

    path2 = mktempdir(cleanup=true)
    model_no_sys =
        DecisionModel(template, sys; optimizer=HiGHS_optimizer, system_to_file=false)

    @test build!(model_no_sys; output_dir=path2) == PSI.BuildStatus.BUILT
    @test solve!(model_no_sys) == RunStatus.SUCCESSFUL

    file_list = sort!(collect(readdir(path2)))
    @test .!all(occursin.(r".h5", file_list))
    ED3 = DecisionModel(path2, HiGHS_optimizer; system=sys)
    build!(ED3, output_dir=path2)
    psi_checksolve_test(ED3, [MOI.OPTIMAL], 240000.0, 10000)
end

@testset "Test NonSpinning reseve model" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc_non_spin", add_reserves=true)
    template = get_thermal_standard_uc_template()
    set_device_model!(
        template,
        DeviceModel(ThermalMultiStart, ThermalStandardUnitCommitment),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserveNonSpinning, NonSpinningReserve, "NonSpinningReserve"),
    )

    UC = DecisionModel(template, c_sys5, optimizer=HiGHS_optimizer)
    output_dir = mktempdir(cleanup=true)
    @test build!(UC; output_dir=output_dir) == PSI.BuildStatus.BUILT
    @test solve!(UC) == RunStatus.SUCCESSFUL
    res = ProblemResults(UC)
    @test isapprox(get_objective_value(res), 247448.0; atol=10000.0)
    vars = res.variable_values
    service_key = PSI.VariableKey(
        ActivePowerReserveVariable,
        PSY.VariableReserveNonSpinning,
        "NonSpinningReserve",
    )
    @test service_key in keys(vars)
end

@testset "Test serialization/deserialization of DecisionModel results" begin
    path = mktempdir(cleanup=true)
    sys = PSB.build_system(PSITestSystems, "c_sys5_re")
    template = get_template_dispatch_with_network(
        NetworkModel(CopperPlatePowerModel; duals=[CopperPlateBalanceConstraint]),
    )
    model = DecisionModel(template, sys; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=path) == PSI.BuildStatus.BUILT
    @test solve!(model, export_problem_results=true) == RunStatus.SUCCESSFUL
    results1 = ProblemResults(model)
    var1_a = read_variable(results1, ActivePowerVariable, ThermalStandard)
    # Ensure that we can deserialize strings into keys.
    var1_b = read_variable(results1, "ActivePowerVariable__ThermalStandard")

    # Results were automatically serialized here.
    results2 = ProblemResults(PSI.get_output_dir(model))
    var2 = read_variable(results2, ActivePowerVariable, ThermalStandard)
    @test var1_a == var2

    # Serialize to a new directory with the exported function.
    results_path = joinpath(path, "results")
    serialize_results(results1, results_path)
    @test isfile(joinpath(results_path, PSI._PROBLEM_RESULTS_FILENAME))
    results3 = ProblemResults(results_path)
    var3 = read_variable(results3, ActivePowerVariable, ThermalStandard)
    @test var1_a == var3
    @test get_system(results3) === nothing
    set_system!(results3, get_system(results1))
    @test get_system(results3) !== nothing

    exp_file =
        joinpath(path, "results", "variables", "ActivePowerVariable__ThermalStandard.csv")
    var4 = PSI.read_dataframe(exp_file)
    # Manually Multiply by the base power var1_a has natural units and export writes directly from the solver
    @test var1_a[:, propertynames(var1_a) .!= :DateTime] == var4 .* 100.0

    @test length(readdir(export_realized_results(results1))) === 6
end

@testset "Test Numerical Stability of Constraints" begin
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    valid_bounds = (coefficient=(min=1.0, max=1.0), rhs=(min=0.4, max=9.930296584))
    model = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

    bounds = PSI.get_constraint_numerical_bounds(model)
    _check_constraint_bounds(bounds, valid_bounds)

    model_bounds = PSI.get_detailed_constraint_numerical_bounds(model)
    valid_model_bounds = Dict(
        :CopperPlateBalanceConstraint__System => (
            coefficient=(min=1.0, max=1.0),
            rhs=(min=6.434489705000001, max=9.930296584),
        ),
        :ActivePowerVariableLimitsConstraint__ThermalStandard__lb =>
            (coefficient=(min=1.0, max=1.0), rhs=(min=Inf, max=-Inf)),
        :ActivePowerVariableLimitsConstraint__ThermalStandard__ub =>
            (coefficient=(min=1.0, max=1.0), rhs=(min=0.4, max=6.0)),
    )
    for (constriant_key, constriant_bounds) in model_bounds
        _check_constraint_bounds(
            constriant_bounds,
            valid_model_bounds[PSI.encode_key(constriant_key)],
        )
    end
end

@testset "Test Numerical Stability of Variables" begin
    template = get_template_basic_uc_simulation()
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc")
    valid_bounds = (min=0.0, max=6.0)
    model = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

    bounds = PSI.get_variable_numerical_bounds(model)
    _check_variable_bounds(bounds, valid_bounds)

    model_bounds = PSI.get_detailed_variable_numerical_bounds(model)
    valid_model_bounds = Dict(
        :StopVariable__ThermalStandard => (min=0.0, max=1.0),
        :StartVariable__ThermalStandard => (min=0.0, max=1.0),
        :ActivePowerVariable__ThermalStandard => (min=0.4, max=6.0),
        :OnVariable__ThermalStandard => (min=0.0, max=1.0),
    )
    for (variable_key, variable_bounds) in model_bounds
        _check_variable_bounds(
            variable_bounds,
            valid_model_bounds[PSI.encode_key(variable_key)],
        )
    end
end

@testset "Decision Model initial_conditions test for ThermalGen" begin
    ######## Test with ThermalStandardUnitCommitment ########
    template = get_thermal_standard_uc_template()
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib"; force_build=true)
    set_device_model!(template, ThermalMultiStart, ThermalStandardUnitCommitment)
    model = DecisionModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    check_duration_on_initial_conditions_values(model, ThermalStandard)
    check_duration_off_initial_conditions_values(model, ThermalStandard)
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with ThermalMultiStartUnitCommitment ########
    template = get_thermal_standard_uc_template()
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib"; force_build=true)
    set_device_model!(template, ThermalMultiStart, ThermalMultiStartUnitCommitment)
    model = DecisionModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT

    check_duration_on_initial_conditions_values(model, ThermalStandard)
    check_duration_off_initial_conditions_values(model, ThermalStandard)
    check_duration_on_initial_conditions_values(model, ThermalMultiStart)
    check_duration_off_initial_conditions_values(model, ThermalMultiStart)
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with ThermalCompactUnitCommitment ########
    template = get_thermal_standard_uc_template()
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib"; force_build=true)
    set_device_model!(template, ThermalMultiStart, ThermalCompactUnitCommitment)
    set_device_model!(template, ThermalStandard, ThermalCompactUnitCommitment)
    model = DecisionModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    check_duration_on_initial_conditions_values(model, ThermalStandard)
    check_duration_off_initial_conditions_values(model, ThermalStandard)
    check_duration_on_initial_conditions_values(model, ThermalMultiStart)
    check_duration_off_initial_conditions_values(model, ThermalMultiStart)
    @test solve!(model) == RunStatus.SUCCESSFUL
end

@testset "Decision Model initial_conditions test for Storage" begin
    ######## Test with BookKeeping ########
    template = get_thermal_dispatch_template_network()
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat"; force_build=true)
    set_device_model!(template, GenericBattery, BookKeeping)
    model = DecisionModel(template, c_sys5_bat; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    check_energy_initial_conditions_values(model, GenericBattery)
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with BatteryAncillaryServices ########
    template = get_thermal_dispatch_template_network()
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat"; force_build=true)
    set_device_model!(template, GenericBattery, BatteryAncillaryServices)
    model = DecisionModel(template, c_sys5_bat; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    check_energy_initial_conditions_values(model, GenericBattery)
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with EnergyTarget ########
    template = get_thermal_dispatch_template_network()
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems"; force_build=true)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    model = DecisionModel(template, c_sys5_bat; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    check_energy_initial_conditions_values(model, BatteryEMS)
    @test solve!(model) == RunStatus.SUCCESSFUL
end

@testset "Decision Model initial_conditions test for Hydro" begin
    ######## Test with HydroDispatchRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd"; force_build=true)
    set_device_model!(template, HydroDispatch, HydroDispatchRunOfRiver)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    model = DecisionModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test !PSI.has_initial_condition_value(
        initial_conditions_data,
        ActivePowerVariable(),
        HydroEnergyReservoir,
    )
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd"; force_build=true)
    set_device_model!(template, HydroDispatch, HydroCommitmentRunOfRiver)
    set_device_model!(template, HydroEnergyReservoir, HydroCommitmentRunOfRiver)
    model = DecisionModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)

    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test PSI.has_initial_condition_value(
        initial_conditions_data,
        OnVariable(),
        HydroEnergyReservoir,
    )
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroDispatchReservoirBudget ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd"; force_build=true)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirBudget)
    model = DecisionModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test !PSI.has_initial_condition_value(
        initial_conditions_data,
        ActivePowerVariable(),
        HydroEnergyReservoir,
    )
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentReservoirBudget ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd"; force_build=true)
    set_device_model!(template, HydroEnergyReservoir, HydroCommitmentReservoirBudget)
    model = DecisionModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test PSI.has_initial_condition_value(
        initial_conditions_data,
        OnVariable(),
        HydroEnergyReservoir,
    )
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroDispatchReservoirStorage ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems"; force_build=true)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    model = DecisionModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test !PSI.has_initial_condition_value(
        initial_conditions_data,
        ActivePowerVariable(),
        HydroEnergyReservoir,
    )
    check_energy_initial_conditions_values(model, HydroEnergyReservoir)
    @test solve!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentReservoirStorage ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems"; force_build=true)
    set_device_model!(template, HydroEnergyReservoir, HydroCommitmentReservoirStorage)
    model = DecisionModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test PSI.has_initial_condition_value(
        initial_conditions_data,
        OnVariable(),
        HydroEnergyReservoir,
    )
    check_energy_initial_conditions_values(model, HydroEnergyReservoir)
    @test solve!(model) == RunStatus.SUCCESSFUL
end

@testset "Test serialization of InitialConditionsData" begin
    sys = PSB.build_system(PSITestSystems, "c_sys5")
    template = get_thermal_standard_uc_template()
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "test"),
    )
    optimizer = GLPK_optimizer

    # Construct and build with default behavior that builds initial conditions.
    model = DecisionModel(template, sys; optimizer=optimizer)
    output_dir = mktempdir(cleanup=true)

    @test build!(model; output_dir=output_dir) == PSI.BuildStatus.BUILT
    ic_file = PSI.get_initial_conditions_file(model)
    test_ic_serialization_outputs(model, ic_file_exists=true, message="make")
    @test solve!(model) == RunStatus.SUCCESSFUL

    # Build again. Initial conditions should be rebuilt.
    PSI.reset!(model)
    @test build!(model; output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=true, message="make")
    @test solve!(model) == RunStatus.SUCCESSFUL

    # Build again, use existing initial conditions.
    model = DecisionModel(
        template,
        sys;
        optimizer=optimizer,
        deserialize_initial_conditions=true,
    )
    @test build!(model; output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=true, message="deserialize")
    @test solve!(model) == RunStatus.SUCCESSFUL

    # Construct and build again with custom initial conditions file.
    initialization_file = joinpath(output_dir, ic_file * ".old")
    mv(ic_file, initialization_file)
    touch(ic_file)
    model = DecisionModel(
        template,
        sys;
        optimizer=optimizer,
        initialization_file=initialization_file,
        deserialize_initial_conditions=true,
    )
    @test build!(model; output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=true, message="deserialize")
    @test solve!(model) == RunStatus.SUCCESSFUL

    # Construct and build again while skipping build of initial conditions.
    rm(ic_file)
    model = DecisionModel(template, sys; optimizer=optimizer, initialize_model=false)
    @test build!(model; output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=false, message="skip")
    @test solve!(model) == RunStatus.SUCCESSFUL

    # Conflicting inputs
    model = DecisionModel(
        template,
        sys;
        optimizer=optimizer,
        initialize_model=false,
        deserialize_initial_conditions=true,
    )
    @test build!(model; output_dir=output_dir, console_level=Logging.AboveMaxLevel) ==
          PSI.BuildStatus.FAILED
    model = DecisionModel(
        template,
        sys;
        optimizer=optimizer,
        initialize_model=false,
        initialization_file="init_file.bin",
    )
    build!(model; output_dir=output_dir, console_level=Logging.AboveMaxLevel) ==
    PSI.BuildStatus.FAILED
end

@testset "Solve with detailed optimizer stats" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    template = get_thermal_standard_uc_template()
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "test"),
    )
    UC = DecisionModel(
        template,
        c_sys5;
        optimizer=GLPK_optimizer,
        detailed_optimizer_stats=true,
    )
    output_dir = mktempdir(cleanup=true)
    @test build!(UC; output_dir=output_dir) == PSI.BuildStatus.BUILT
    @test solve!(UC) == RunStatus.SUCCESSFUL
    # We only test this field because most free solvers don't support detailed stats
    @test !ismissing(get_optimizer_stats(UC).objective_bound)
end


text:
using TimeSeriesClustering
using JLD2
using Cbc
using Random

reference_results = Dict{String,Any}()

Random.seed!(1111)
for data in [:CEP_GER1,:CEP_GER18]
    ts_input_data = load_timeseries_data(data)
    #mr: method, representation, n_init
    mr = [["kmeans","centroid",1000],
    ["kmeans","medoid",1000],
    ["kmedoids","centroid",1000],
    ["kmedoids","medoid",1000],
    ["hierarchical","centroid",1],
    ["hierarchical","medoid",1]]
    # default
    for (method,repr,n_init) in mr
        Random.seed!(1111)
        try
            reference_results["$data-$method-$repr-default"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=n_init)
        catch
            reference_results["$data-$method-$repr-default"] = "not defined"
            println("$data-$method-$repr-default not defined")
        end
    end
    # n_clust=1
    for (method,repr,n_init) in mr
        Random.seed!(1111)
        try
            reference_results["$data-$method-$repr-n_clust1"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=1,n_init=n_init)
        catch
            reference_results["$data-$method-$repr-n_clust1"] = "not defined"
            println("$data-$method-$repr-n_clust1 not defined")
        end
    end
    # n_clust = N
    for (method,repr,n_init) in mr
        Random.seed!(1111)
        try
            reference_results["$data-$method-$repr-n_clustK"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=ts_input_data.K,n_init=n_init)
        catch
            reference_results["$data-$method-$repr-n_clustK"] = "not defined"
            println("$data-$method-$repr-n_clustK not defined")
        end
    end
end

data = :CEP_GER1
ts_input_data = load_timeseries_data(data)
method = "kmedoids_exact"
repr = "medoid"
using Cbc
optimizer = Cbc.Optimizer
reference_results["$data-$method-$repr-default"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=1,kmexact_optimizer=optimizer)
reference_results["$data-$method-$repr-n_clust1"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=1,n_init=1,kmexact_optimizer=optimizer)
reference_results["$data-$method-$repr-n_clustK"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=ts_input_data.K,n_init=1,kmexact_optimizer=optimizer)

# MultiClustAtOnce
Random.seed!(1111)
method = "hierarchical"
repr = "centroid"
reference_results["$data-$method-$repr-MultiClust"] = run_clust(ts_input_data,[1,5,ts_input_data.K];method=method,representation=repr,n_init=1)

# ClustResultAll
Random.seed!(1111)
method = "hierarchical"
repr = "medoid"
reference_results["$data-$method-$repr-ClustResultAll"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=10,get_all_clust_results=true)

# attribute weighting
method = "hierarchical"
repr = "centroid"
attribute_weights=Dict("solar"=>1.0, "wind"=>2.0, "el_demand"=>3.0)
reference_results["$data-$method-$repr-AttributeWeighting"] = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=1,attribute_weights=attribute_weights)

@save normpath(joinpath(dirname(@__FILE__),"run_clust.jld2")) reference_results

ext: jl


text:
module MyModule
export hello

hello() = "Hello my module"

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 28761bd0-f1fb-11ea-1f82-717a48c220a1
using Pkg

# ╔═╡ 6292b148-f1fb-11ea-2c32-91cd70cb3bd9
function useit(list::Array{Symbol})
    installed = [key for key in keys(Pkg.installed())]
    strpackages = @. string(list)
    uninstalled = setdiff(strpackages,installed)

    map(Pkg.add,uninstalled)
    for package ∈ list
        @eval using $package
    end
end

# ╔═╡ 036d9cfe-f1fc-11ea-2b6b-4fae4cf7dc02
Pkg.add(url="https://github.com/loopulos/NetABM.jl")

# ╔═╡ 62c3033e-f1fb-11ea-32e1-e94bbd4058f6
thep2 = [:LightGraphs, :StatsBase, :Plots, :GraphPlot, :Statistics, :LaTeXStrings, :Bootstrap]

# ╔═╡ 62c3d066-f1fb-11ea-3320-abbb1d14df89
useit(thep2)

# ╔═╡ 5eaa7916-f1fc-11ea-3354-fb5b33917bc7
global rtall = Array{Array}(undef,0);

# ╔═╡ 75ae33b4-f1fc-11ea-0ddb-9f912fed94e9
global lrtall = Array{Array}(undef,0);

# ╔═╡ 75b0a64e-f1fc-11ea-3245-ade5236eabf4
global raall = Array{Array}(undef,0);

# ╔═╡ 75c10a8e-f1fc-11ea-2949-5784fe7301d2
global lraall = Array{Array}(undef,0);

# ╔═╡ 75ebbe1e-f1fc-11ea-2d02-9fed550f2214
global n_boot = 100

# ╔═╡ 76001a58-f1fc-11ea-277c-e9e68b50fb9a
global cil = 0.95

# ╔═╡ 7600de2a-f1fc-11ea-07d8-c32270c2e0bd
global j

# ╔═╡ 76353418-f1fc-11ea-0fa4-73369a86b186
@time let
    #  for j in [0.0, 0.1, 0.4, 0.7, 1.0]
    for j in [0.5]
        Imean = Array{Float64}(undef,0)
        Ilow = Array{Float64}(undef,0)
        Ihigh = Array{Float64}(undef,0)
        Iall = Array{Array}(undef,0)
        for a in 1:100
            g = barabasi_albert(num_agents,3)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            #  get_coop!(agents);
            I = Array{Float64}(undef,0);
            rt = Array{Float64}(undef,0);
            ra = Array{Float64}(undef,0);
            lrt = Array{Float64}(undef,0);
            lra = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=false);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                #  get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
                push!(rt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && !ag.adapter]));
                push!(lrt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && ag.adapter]));
                push!(ra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && !ag.adapter]));
                push!(lra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && ag.adapter]));
            end
            push!(Iall,I)
            push!(rtall,rt)
            push!(lrtall,lrt)
            push!(raall,ra)
            push!(lraall,lra)
            GC.gc()
        end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
            #  GC.gc()
        end
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
    GC.gc()
end

# ╔═╡ 764b35e2-f1fc-11ea-0d88-630a537082e9
Imean = Array{Float64}(undef,0)

# ╔═╡ 764c297a-f1fc-11ea-110c-8181fb96b08c
Ilow = Array{Float64}(undef,0)

# ╔═╡ 76686946-f1fc-11ea-0fa2-b7edcca389b4
Ihigh = Array{Float64}(undef,0)

# ╔═╡ 7682b67a-f1fc-11ea-290e-3be1c6b72b2c
Iall = Array{Array}(undef,0)

# ╔═╡ 76840afc-f1fc-11ea-33d8-e7a199ac97f6
for step in eachindex(lraall[1])
    bs1 = bootstrap(mean,getindex.(lraall,step),BasicSampling(n_boot))
    bs2 = bootstrap(std,getindex.(lraall,step),BasicSampling(n_boot))
    bci1 = confint(bs1, BasicConfInt(cil))
    bci2 = confint(bs2, BasicConfInt(cil))
    push!(Imean,bci1[1][1])
    push!(Ilow,bci2[1][2])
    push!(Ihigh,bci2[1][3])
end

# ╔═╡ 76a007c0-f1fc-11ea-04cc-0d91a50209ff
pyplot()

# ╔═╡ 76cb3094-f1fc-11ea-347d-c74f8474d987
p

# ╔═╡ 8642ec88-f1fc-11ea-0b4c-1bea3b95596e
begin
		num_agents = 1000
		#  g = erdos_renyi(num_agents,0.4)
		g = barabasi_albert(num_agents,9)
		agents = [Agent(i) for i = 1:num_agents]
		init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
		set_coop_agents!(agents;p_cop=0.1)
		map(x -> assign_contacts!(g,x), agents);
		get_coop!(agents)
end

# ╔═╡ 76041c34-f1fc-11ea-16a6-3987f38e1d5b
global adoptstep = (0,[0.1])

# ╔═╡ 761bee22-f1fc-11ea-0a83-31849f4d3661
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")

# ╔═╡ 76b68ab8-f1fc-11ea-2aaa-711f0ef72d4f
p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label="lrt")

# ╔═╡ 75d8b670-f1fc-11ea-3505-996f9c062d6c
global num_agents = 1000

# ╔═╡ 761ab8f4-f1fc-11ea-303b-4d26188526b2
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")

# ╔═╡ 75c18024-f1fc-11ea-2fcd-a91bcc345fc0
global adoptstep = (1,[0.1,0.05])

# ╔═╡ Cell order:
# ╠═28761bd0-f1fb-11ea-1f82-717a48c220a1
# ╠═036d9cfe-f1fc-11ea-2b6b-4fae4cf7dc02
# ╠═6292b148-f1fb-11ea-2c32-91cd70cb3bd9
# ╠═62c3033e-f1fb-11ea-32e1-e94bbd4058f6
# ╠═62c3d066-f1fb-11ea-3320-abbb1d14df89
# ╠═8642ec88-f1fc-11ea-0b4c-1bea3b95596e
# ╠═5eaa7916-f1fc-11ea-3354-fb5b33917bc7
# ╠═75ae33b4-f1fc-11ea-0ddb-9f912fed94e9
# ╠═75b0a64e-f1fc-11ea-3245-ade5236eabf4
# ╠═75c10a8e-f1fc-11ea-2949-5784fe7301d2
# ╠═75c18024-f1fc-11ea-2fcd-a91bcc345fc0
# ╠═75d8b670-f1fc-11ea-3505-996f9c062d6c
# ╠═75ebbe1e-f1fc-11ea-2d02-9fed550f2214
# ╠═76001a58-f1fc-11ea-277c-e9e68b50fb9a
# ╠═7600de2a-f1fc-11ea-07d8-c32270c2e0bd
# ╠═76041c34-f1fc-11ea-16a6-3987f38e1d5b
# ╠═761ab8f4-f1fc-11ea-303b-4d26188526b2
# ╠═761bee22-f1fc-11ea-0a83-31849f4d3661
# ╠═76353418-f1fc-11ea-0fa4-73369a86b186
# ╠═764b35e2-f1fc-11ea-0d88-630a537082e9
# ╠═764c297a-f1fc-11ea-110c-8181fb96b08c
# ╠═76686946-f1fc-11ea-0fa2-b7edcca389b4
# ╠═7682b67a-f1fc-11ea-290e-3be1c6b72b2c
# ╠═76840afc-f1fc-11ea-33d8-e7a199ac97f6
# ╠═76a007c0-f1fc-11ea-04cc-0d91a50209ff
# ╠═76b68ab8-f1fc-11ea-2aaa-711f0ef72d4f
# ╠═76cb3094-f1fc-11ea-347d-c74f8474d987


text:
struct TestLasPoint
    x::Int32
    y::Int32
    z::Int32
    intensity::UInt16
    flag_byte::UInt8  # return number (3 bits), number of returns (3 bits), scan direction flag, edge of flight line
    raw_classification::UInt8  # classification (5 bits), classification flags
    flag_byte_1::UInt8  # return number (4 bits) & number of returns (4 bits)
    flag_byte_2::UInt8 # classification flags, scanner channel,
    classification::UInt8
    user_data::UInt8
    scan_angle::Int8
    extended_scan_angle::Int16
    pt_src_id::UInt16
    gps_time::Float64
    red::N0f16
    green::N0f16
    blue::N0f16
    nir::N0f16
    wave_packet_descriptor_index::UInt8
    wave_packet_byte_offset::UInt64
    wave_packet_size_in_bytes::UInt32
    wave_return_location::Float32
    wave_x_t::Float32
    wave_y_t::Float32
    wave_z_t::Float32
end

function TestLasPoint()
    return TestLasPoint(
        rand(Int32), # x
        rand(Int32), # y
        rand(Int32), # z
        rand(UInt16), # intensity::UInt16
        rand(UInt8), # flag_byte::UInt8  # return number (3 bits), number of returns (3 bits), scan direction flag, edge of flight line
        rand(UInt8), # raw_classification::UInt8  # classification (5 bits), classification flags
        rand(UInt8), # flag_byte_1::UInt8  # return number (4 bits) & number of returns (4 bits)
        rand(UInt8), # flag_byte_2::UInt8 # classification flags, scanner channel,
        rand(UInt8), # classification::UInt8
        rand(UInt8), # user_data::UInt8
        rand(Int8), # scan_angle::Int16
        rand(Int16), # scan_angle::Int16
        rand(UInt16), # pt_src_id::UInt16
        rand(Float64), # gps_time::Float64
        rand(N0f16), # red::N0f16
        rand(N0f16), # green::N0f16
        rand(N0f16), # blue::N0f16
        rand(N0f16), # nir::N0f16
        rand(UInt8), # wave_packet_descriptor_index::UInt8
        rand(UInt64), # wave_packet_byte_offset::UInt64
        rand(UInt32), # wave_packet_size_in_bytes::UInt32
        rand(Float32), # wave_return_location::Float32
        rand(Float32), # wave_x_t::Float32
        rand(Float32), # wave_y_t::Float32
        rand(Float32), # wave_z_t::Float32
    )
end

function Base.:(==)(h1::T, h2::T) where T <:Union{LasHeader, LasPoint}
    if typeof(h1) != typeof(h2)
        return false
    end
    header_fields = fieldnames(typeof(h1))
    if all([getproperty(h1, p) == getproperty(h2, p) for p in header_fields])
        return true
    else
        @error "The following LasHeader values did not match: $(header_fields[[getproperty(h1, p) != getproperty(h2, p) for p in header_fields]])"
        return false
    end
end

# point formats 0 - 5
get_constructor_values(::Type{LasPoint0}, t::TestLasPoint) = [t.x, t.y, t.z, t.intensity, t.flag_byte, t.raw_classification, t.scan_angle, t.user_data, t.pt_src_id]
get_constructor_values(::Type{LasPoint1}, t::TestLasPoint) = [get_constructor_values(LasPoint0, t)..., t.gps_time]
get_constructor_values(::Type{LasPoint2}, t::TestLasPoint) = [get_constructor_values(LasPoint0, t)..., t.red, t.green, t.blue]
get_constructor_values(::Type{LasPoint3}, t::TestLasPoint) = [get_constructor_values(LasPoint0, t)..., t.gps_time, t.red, t.green, t.blue]
get_constructor_values(::Type{LasPoint4}, t::TestLasPoint) = [get_constructor_values(LasPoint1, t)..., t.wave_packet_descriptor_index, t.wave_packet_byte_offset, t.wave_packet_size_in_bytes, t.wave_return_location, t.wave_x_t, t.wave_y_t, t.wave_z_t]
get_constructor_values(::Type{LasPoint5}, t::TestLasPoint) = [get_constructor_values(LasPoint3, t)..., t.wave_packet_descriptor_index, t.wave_packet_byte_offset, t.wave_packet_size_in_bytes, t.wave_return_location, t.wave_x_t, t.wave_y_t, t.wave_z_t]

# point formats 6 - 10
get_constructor_values(::Type{LasPoint6}, t::TestLasPoint) = [t.x, t.y, t.z, t.intensity, t.flag_byte_1, t.flag_byte_2, t.classification, t.user_data, t.extended_scan_angle, t.pt_src_id, t.gps_time]
get_constructor_values(::Type{LasPoint7}, t::TestLasPoint) = [get_constructor_values(LasPoint6, t)..., t.red, t.green, t.blue]
get_constructor_values(::Type{LasPoint8}, t::TestLasPoint) = [get_constructor_values(LasPoint7, t)..., t.nir]
get_constructor_values(::Type{LasPoint9}, t::TestLasPoint) = [get_constructor_values(LasPoint6, t)..., t.wave_packet_descriptor_index, t.wave_packet_byte_offset, t.wave_packet_size_in_bytes, t.wave_return_location, t.wave_x_t, t.wave_y_t, t.wave_z_t]
get_constructor_values(::Type{LasPoint10}, t::TestLasPoint) = [get_constructor_values(LasPoint8, t)..., t.wave_packet_descriptor_index, t.wave_packet_byte_offset, t.wave_packet_size_in_bytes, t.wave_return_location, t.wave_x_t, t.wave_y_t, t.wave_z_t]

const SUPPORTED_LAS_POINT_FORMATS = [
    LasPoint0,
    LasPoint1,
    LasPoint2,
    LasPoint3,
    LasPoint4, # with waveform
    LasPoint5, # with waveform
    LasPoint6,
    LasPoint7,
    LasPoint8,
    LasPoint9, # with waveform
    LasPoint10, # with waveform
]

@testset "Test different LAS Versions" begin
    workdir = mktempdir()
    for version_minor = 1 : 4
        for point_type_index = 1 : length(SUPPORTED_LAS_POINT_FORMATS)
            point_type = SUPPORTED_LAS_POINT_FORMATS[point_type_index]
            number_of_points = 10
            max_return_number = version_minor > 3 ? 15 : 5
            return_numbers = rand(1 : max_return_number, max_return_number)
            return_count = [UInt32(sum(i .== return_numbers)) for i = 1 : 5]

            # use LasPoint1
            points = [point_type(get_constructor_values(point_type, TestLasPoint())...) for i = 1 : number_of_points]
            data_record_length = size(point_type)

            point_return_count = version_minor == 4 ? zeros(UInt64, 15) : zeros(UInt32, 5)
            return_number_array = clamp.([return_number(p) for p in points], 1, length(point_return_count))
            @inbounds for i in 1:length(point_return_count)
                point_return_count[i] = sum(return_number_array .== i)
            end

            scale = [0.001, 0.002, 0.003]
            offset = [1.0, 2.0, 3.0]
            x_min, x_max = extrema([p.x * scale[1] + offset[1] for  p in points])
            y_min, y_max = extrema([p.y * scale[2] + offset[2] for  p in points])
            z_min, z_max = extrema([p.z * scale[3] + offset[3] for  p in points])

            header_size = 227
            if version_minor == 3
                header_size = 235
            elseif version_minor == 4
                header_size = 375
            end
            t = Dates.now()
            header = LasHeader(
                UInt16(0), # file_source_id
                UInt16(0), # global_encoding
                UInt32(0), # guid_1
                UInt16(0), # guid_2
                UInt16(0), # guid_3
                "", # guid_4
                UInt8(1), # LAS version major
                UInt8(version_minor), # LAS version minor
                "OTHER", # System
                "Julia LasIO writer", # Software
                UInt16(dayofyear(t)), # creation_year
                UInt16(year(t)), # creation_dayofyear
                version_minor == 4 ? UInt16(header_size) : UInt16(header_size), # header_size
                version_minor == 4 ? UInt32(header_size) : UInt32(header_size), # data_offset
                UInt32(0), # n_vlr
                UInt8(point_type_index-1), # data_format_id
                UInt16(data_record_length), # data_record_length
                version_minor == 4 ? UInt32(0) : number_of_points,
                version_minor == 4 ? fill(UInt32(0), 5) : point_return_count,
                0.001,
                0.001,
                0.001,
                0.0,
                0.0,
                0.0,
                x_max,
                x_min,
                y_max,
                y_min,
                z_max,
                z_min,
                UInt64(0), # start_of_waveform_data_packet_record
                UInt64(0), # start_of_first_extended_variable_length_record
                UInt32(0), # number_of_extended_variable_length_records
                version_minor == 4 ? number_of_points : UInt64(0),  # extended_number_of_point_records
                version_minor == 4 ? point_return_count : fill(UInt32(0), 15), # extended_number_of_points_by_return
                [], # variable_length_records
                [], # user_defined_bytes
            )

            # save data to LAS file
            las_output_file_path = joinpath(workdir, "test_version_1_$(version_minor)_pointformat_LasPoint$(point_type_index-1).las")
            save(las_output_file_path, header, points)

            # read las file
            las_input_header, las_points_read = load(las_output_file_path)
            @test las_input_header == header
            @test all(points .== las_points_read)

            # save data to LAZ file
            laz_output_file_path = joinpath(workdir, "test_version_1_$(version_minor)_pointformat_LasPoint$(point_type_index-1).laz")
            save(laz_output_file_path, header, points)

            # read laz file
            laz_input_header, laz_points_read = load(las_output_file_path)
            @test laz_input_header == header
            @test all(points .== laz_points_read)

        end
    end
    rm(workdir, recursive=true)
end


text:
function create(c::Client, ::Type{Reaction}, m::Message, emoji::StringOrChar)
    return create_reaction(c, m.channel_id, m.id, emoji)
end
function create(c::Client, ::Type{Reaction}, m::Message, e::Emoji)
    emoji = e.id === nothing ? e.name : "$(e.name):$(e.id)"
    return create_reaction(c, m.channel_id, m.id, emoji)
end

function retrieve(c::Client, ::Type{Reaction}, m::Message, emoji::StringOrChar)
    return get_reactions(c, m.channel_id, m.id, emoji)
end
retrieve(c::Client, ::Type{Reaction}, m::Message, e::Emoji) = get(Reaction, c, m, e.name)

function delete(c::Client, ::Type{Reaction}, m::Message)
    return delete_all_reactions(c, m.channel_id, m.id)
end
function delete(c::Client, emoji::StringOrChar, m::Message, u::User)
    return if c.state.user !== nothing && u.id == c.state.user.id
        delete_own_reaction(c, m.channel_id, m.id, emoji)
    else
        delete_user_reaction(c, m.channel_id, m.id, emoji, u.id)
    end
end
function delete(c::Client, e::Emoji, m::Message, u::User)
    emoji = e.id === nothing ? e.name : "$(e.name):$(e.id)"
    delete(c, emoji, m, u)
end


text: struct NeptuneType <: abstractPlanet end


text:
module ZeroRPC

import ZMQ
import Msgpack
import Base.Random

function unpack(envelope::Array{ZMQ.Message})
  raw_msg = convert(IOStream, envelope[3])
  seek(raw_msg, 0)
  return Msgpack.unpack(raw_msg.data)
end

function generate_heartbeat(response_to)
  header = { "v" => 3, "message_id" => string(Random.uuid4()), "response_to" => response_to}
  response = Any[header, "_zpc_hb", [0]]
  raw_response = ZMQ.Message(Msgpack.pack(response))
  return raw_response
end

type Server
  functions::Dict{String, Function}
  ctx::ZMQ.Context
  socket::ZMQ.Socket

  function Server(_module::Module)

    functions = Dict{String, Function}()

    for func_symbol in names(_module)
      f = eval(_module, func_symbol)
      if isa(f, Function)
        functions[string(func_symbol)] = f
      end
    end

    ctx = ZMQ.Context(1)
    socket = ZMQ.Socket(ctx, ZMQ.ROUTER)
    return new (functions, ctx, socket)
  end
end

function bind(server::Server, dest)
  ZMQ.bind(server.socket, dest)
end

function run(server::Server)
  while true
    envelope::Array{ZMQ.Message} = ZMQ.recv_multipart(server.socket)
    message = unpack(envelope)

    if message[2] == "_zpc_hb"
      # ignore incoming heartbeats for now
      continue
    end

    response_type = "OK"
    ret = nothing
    f = function()
          try
            ret = Any[apply(server.functions[message[2]], message[3])]
          catch e
            buf = IOBuffer()
            showerror(buf, e, catch_backtrace())
            bt = takebuf_string(buf)
            println(bt)
            ret = Any[string(e), string(e), bt]
            response_type = "ERR"
          end
        end
    task = @async f()
    a = time()
    while !istaskdone(task)
      b = time()
      if b - a > 5
        # send heartbeat
        env = deepcopy(envelope)
        env[3] = generate_heartbeat(message[1]["message_id"])
        ZMQ.send_multipart(server.socket, env)
        a = b
      end
      yield()
    end

    header = { "v" => 3, "message_id" => string(Random.uuid4()), "response_to" => message[1]["message_id"] }
    response = Any[header, response_type, ret]
    raw_response = ZMQ.Message(Msgpack.pack(response))
    envelope[3] = raw_response
    ZMQ.send_multipart(server.socket, envelope)
  end
end

end # module


text:
server="ai.ku.edu.tr/"
if !isdir("data/demo")
    println("Downloading sample questions and images from CLEVR dataset...")
    download(server*"data/mac-network/demo.tar.gz","demo.tar.gz")
    run(`tar -xzf demo.tar.gz`)
    rm("demo.tar.gz")
end

if !isfile("models/macnet2.jld2")
    println("Downloading pre-trained model from our servers...")
    download(server*"models/mac-network/macnet2.jld2","models/macnet2.jld2")
end
println("Demo setup is completed")


text:
include("helmholtzhodge.jl")

using Random
Random.seed!(0)

function sphrandn(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, 2n-1)
    for i = 1:m
        A[i,1] = randn(T)
    end
    for j = 1:n-1
        for i = 1:m-j
            A[i,2j] = randn(T)
            A[i,2j+1] = randn(T)
        end
    end
    A
end

#=
N = 2 .^(5:13)
t = zeros(length(N), 2)
err = zeros(length(N))

j = 1
for n in N
    println("This is n: ", n)
    U1 = sphrandn(Float64, n, n)
    U2 = sphrandn(Float64, n, n)
    U1[1] = 0
    U2[1] = 0
    Us = zero(U1); Ut = zero(U1);
    V1 = zero(U1); V2 = zero(U1); V3 = zero(U1); V4 = zero(U1);

    HH = HelmholtzHodge(Float64, n)
    t[j, 1] = @elapsed for k = 1:10
        HelmholtzHodge(Float64, n)
    end
    t[j, 1] /= 10

    gradient!(U1, V1, V2)
    curl!(U2, V3, V4)

    V5 = V1 + V3
    V6 = V2 + V4

    helmholtzhodge!(HH, Us, Ut, V5, V6)

    t[j, 2] = @elapsed for k = 1:10
        helmholtzhodge!(HH, Us, Ut, V5, V6)
    end
    t[j, 2] /= 10
    global j += 1
end

j = 1
for n in N
    println("This is n: ", n)
    HH = HelmholtzHodge(Float64, n)
    for k = 1:10
        U1 = sphrandn(Float64, n, n)
        U2 = sphrandn(Float64, n, n)
        U1[1] = 0
        U2[1] = 0
        Us = zero(U1); Ut = zero(U1);
        V1 = zero(U1); V2 = zero(U1); V3 = zero(U1); V4 = zero(U1);

        gradient!(U1, V1, V2)
        curl!(U2, V3, V4)

        V5 = V1 + V3
        V6 = V2 + V4

        helmholtzhodge!(HH, Us, Ut, V5, V6)

        println("This is the spheroidal relative backward error: ", norm(Us-U1)/norm(U1))
        println("This is the toroidal relative backward error: ", norm(Ut-U2)/norm(U2))
        err[j] += norm(Us-U1)/norm(U1) + norm(Ut-U2)/norm(U2)
    end
    err[j] /= 10
    global j += 1
end
=#

using PyPlot

cd("/Users/Mikael/Dropbox/Helmholtz")

t = [0.000301073 8.59417e-5; 0.000870661 0.000313344; 0.00371459 0.0013411; 0.0158751 0.00519203; 0.0629211 0.0214159; 0.272275 0.0841661; 1.18461 0.338013; 5.22898 1.34888; 20.2509 5.41399]
err = [9.73175783254202e-16, 1.1848048133640504e-15, 1.2864591012075482e-15, 1.4887354817191668e-15, 2.0250488931546365e-15, 2.714940169354596e-15, 3.426270329252532e-15, 4.804155288069517e-15, 6.620111175740055e-15]

clf()
loglog(N.-1, t[:, 1], "xk", N.-1, t[:, 2], "+k", N.-1, (N.-1).^2/2e6, "-k")
xlabel("Degree \$n\$"); ylabel("Execution Time (s)"); grid()
legend(["Pre-computation","Execution","\$\\mathcal{O}(n^2)\$"])
savefig("helmholtzhodgetime.pdf")

clf()
loglog(N.-1, err, "xk", N.-1, sqrt.(N.*log.(N))*eps()/7, "-k")
xlabel("Degree \$n\$"); ylabel("Relative Error"); grid()
ylim((6e-16,1.2e-14))
legend(["Error","\$\\mathcal{O}(\\sqrt{n \\log(n)}\\epsilon)\$"])
savefig("helmholtzhodgeerror.pdf")