text:
function get_deserialized(sim::Simulation, stage_info)
    path = mktempdir()
    directory = PSI.serialize(sim; path = path)
    return Simulation(directory, stage_info)
end

function test_load_simulation(file_path::String)

    c_sys5_uc = build_system("c_sys5_uc")
    single_stage_definition =
        Dict("ED" => Stage(GenericOpProblem, template_ed, c_sys5_uc, ipopt_optimizer))

    single_sequence = SimulationSequence(
        step_resolution = Hour(1),
        order = Dict(1 => "ED"),
        horizons = Dict("ED" => 12),
        intervals = Dict("ED" => (Hour(1), Consecutive())),
        ini_cond_chronology = IntraStageChronology(),
    )

    sim_single = Simulation(
        name = "consecutive",
        steps = 2,
        stages = single_stage_definition,
        stages_sequence = single_sequence,
        simulation_folder = file_path,
    )
    build!(sim_single)
    res = execute!(sim_single)

    @testset "Single stage sequential tests" begin
        stage_single = PSI.get_stage(sim_single, "ED")
        @test JuMP.termination_status(stage_single.internal.psi_container.JuMPmodel) in
              [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    end

    stage_info = Dict(
        "UC" => Dict("optimizer" => GLPK_optimizer, "jump_model" => nothing),
        "ED" => Dict("optimizer" => ipopt_optimizer),
    )
    # Tests of a Simulation without Caches
    duals = [:CopperPlateBalance]
    c_sys5_hy_uc = build_system("c_sys5_hy_uc")
    c_sys5_hy_ed = build_system("c_sys5_hy_ed")
    stages_definition = Dict(
        "UC" => Stage(
            GenericOpProblem,
            template_hydro_basic_uc,
            c_sys5_hy_uc,
            stage_info["UC"]["optimizer"];
            constraint_duals = duals,
        ),
        "ED" => Stage(
            GenericOpProblem,
            template_hydro_ed,
            c_sys5_hy_ed,
            stage_info["ED"]["optimizer"];
            constraint_duals = duals,
        ),
    )

    sequence = SimulationSequence(
        step_resolution = Hour(24),
        order = Dict(1 => "UC", 2 => "ED"),
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        horizons = Dict("UC" => 24, "ED" => 12),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_from_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
            ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                variable_from_stage = PSI.ACTIVE_POWER,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterStageChronology(),
    )
    sim = Simulation(
        name = "aggregation",
        steps = 2,
        stages = stages_definition,
        stages_sequence = sequence,
        simulation_folder = file_path,
    )

    # Run twice, once building normally, once after deserializing.
    for i in 1:2
        output_dir = "test" * string(i)
        if i == 2
            sim = get_deserialized(sim, stage_info)
        end

        build!(sim; output_dir = output_dir, recorders = [:simulation])
        sim_results = execute!(sim)

        stage_names = keys(sim.stages)
        step = ["step-1", "step-2"]

        @testset "All stages executed - No Cache" begin
            for name in stage_names
                stage = PSI.get_stage(sim, name)
                @test JuMP.termination_status(stage.internal.psi_container.JuMPmodel) in
                      [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
            end
        end

        @testset "Test reading and writing to the results folder" begin
            for name in stage_names
                files = collect(readdir(sim_results.results_folder))
                for f in files
                    rm("$(sim_results.results_folder)/$f")
                end
                rm(sim_results.results_folder)
                res = load_simulation_results(sim_results, name)
                !ispath(res.results_folder) && mkdir(res.results_folder)
                write_results(res)
                loaded_res = load_operation_results(sim_results.results_folder)
                @test loaded_res.variable_values == res.variable_values
                @test loaded_res.parameter_values == res.parameter_values
            end
        end

        @testset "Test file names" begin
            for name in stage_names
                files = collect(readdir(sim_results.results_folder))
                for f in files
                    rm("$(sim_results.results_folder)/$f")
                end
                res = load_simulation_results(sim_results, name)
                write_results(res)
                variable_list = String.(PSI.get_variable_names(sim, name))
                variable_list = [
                    variable_list
                    "dual_CopperPlateBalance"
                    "optimizer_log"
                    "time_stamp"
                    "check"
                    "base_power"
                    "parameter_P_InterruptibleLoad"
                    "parameter_P_PowerLoad"
                    "parameter_P_RenewableDispatch"
                    "parameter_P_HydroEnergyReservoir"
                ]
                file_list = collect(readdir(sim_results.results_folder))
                for name in file_list
                    variable = splitext(name)[1]
                    @test any(x -> x == variable, variable_list)
                end
            end
        end

        @testset "Test argument errors" begin
            for name in stage_names
                res = load_simulation_results(sim_results, name)
                if isdir(res.results_folder)
                    files = collect(readdir(res.results_folder))
                    for f in files
                        rm("$(res.results_folder)/$f")
                    end
                    rm("$(res.results_folder)")
                end
                @test_throws IS.ConflictingInputsError write_results(res)
            end
        end

        @testset "Test simulation output serialization and deserialization" begin
            output_path = joinpath(dirname(sim_results.results_folder), "output_references")
            sim_output = collect(readdir(output_path))
            @test sim_output == [
                "base_power.json",
                "chronologies.json",
                "results_folder.json",
                "stage-ED",
                "stage-UC",
            ]
            sim_test = PSI.deserialize_sim_output(dirname(output_path))
            @test sim_test.ref == sim_results.ref
        end

        @testset "Test load simulation results between the two methods of load simulation" begin
            for name in stage_names
                variable = PSI.get_variable_names(sim, name)
                results = load_simulation_results(sim_results, name)
                res = load_simulation_results(sim_results, name, step, variable)
                @test results.variable_values == res.variable_values
                @test results.parameter_values == res.parameter_values
            end
        end

        @testset "Test verify length of time_stamp" begin
            for name in keys(sim.stages)
                results = load_simulation_results(sim_results, name)
                @test size(unique(results.time_stamp), 1) == size(results.time_stamp, 1)
            end
        end

        @testset "Test verify no gaps in the time_stamp" begin
            for name in keys(sim.stages)
                stage = sim.stages[name]
                results = load_simulation_results(sim_results, name)
                resolution = convert(
                    Dates.Millisecond,
                    PSY.get_forecasts_resolution(PSI.get_sys(stage)),
                )
                time_stamp = results.time_stamp
                length = size(time_stamp, 1)
                test = results.time_stamp[1, 1]:resolution:results.time_stamp[length, 1]
                @test time_stamp[!, :Range] == test
            end
        end
        ###########################################################

        @testset "Test dual constraints in results" begin
            res = PSI.load_simulation_results(sim_results, "ED")
            dual =
                JuMP.dual(sim.stages["ED"].internal.psi_container.constraints[:CopperPlateBalance][1])
            @test isapprox(
                dual,
                res.dual_values[:dual_CopperPlateBalance][1, 1],
                atol = 1.0e-4,
            )
            !ispath(res.results_folder) && mkdir(res.results_folder)
            PSI.write_to_CSV(res)
            @test !isempty(res.results_folder)
        end

        @testset "Test verify parameter feedforward for consecutive UC to ED" begin
            P_keys = [
                (PSI.ACTIVE_POWER, PSY.HydroEnergyReservoir),
                #(PSI.ON, PSY.ThermalStandard),
                #(PSI.ACTIVE_POWER, PSY.HydroEnergyReservoir),
            ]

            vars_names = [
                PSI.variable_name(PSI.ACTIVE_POWER, PSY.HydroEnergyReservoir),
                #PSI.variable_name(PSI.ON, PSY.ThermalStandard),
                #PSI.variable_name(PSI.ACTIVE_POWER, PSY.HydroEnergyReservoir),
            ]
            for (ik, key) in enumerate(P_keys)
                variable_ref = PSI.get_reference(sim_results, "UC", 1, vars_names[ik])[1] # 1 is first step
                array = PSI.get_parameter_array(PSI.get_parameter_container(
                    sim.stages["ED"].internal.psi_container,
                    Symbol(key[1]),
                    key[2],
                ))
                parameter = collect(values(value.(array.data)))  # [device, time] 1 is first execution
                raw_result = Feather.read(variable_ref)
                for j in 1:size(parameter, 1)
                    result = raw_result[end, j] # end is last result [time, device]
                    initial = parameter[1] # [device, time]
                    @test isapprox(initial, result)
                end
            end
        end

        @testset "Test verify time gap for Consecutive" begin
            names = ["UC"]
            for name in names
                variable_list = PSI.get_variable_names(sim, name)
                reference_1 = PSI.get_reference(sim_results, name, 1, variable_list[1])[1]
                reference_2 = PSI.get_reference(sim_results, name, 2, variable_list[1])[1]
                time_file_path_1 = joinpath(dirname(reference_1), "time_stamp.feather") #first line, file path
                time_file_path_2 = joinpath(dirname(reference_2), "time_stamp.feather")
                time_1 = convert(Dates.DateTime, Feather.read(time_file_path_1)[end, 1]) # first time
                time_2 = convert(Dates.DateTime, Feather.read(time_file_path_2)[1, 1])
                @test time_2 == time_1
            end
        end

        @testset "Test verify initial condition feedforward for consecutive ED to UC" begin
            ic_keys = [PSI.ICKey(PSI.DevicePower, PSY.ThermalStandard)]
            vars_names = [PSI.variable_name(PSI.ACTIVE_POWER, PSY.ThermalStandard)]
            for (ik, key) in enumerate(ic_keys)
                variable_ref = PSI.get_reference(sim_results, "ED", 1, vars_names[ik])[24]
                initial_conditions =
                    get_initial_conditions(PSI.get_psi_container(sim, "UC"), key)
                for ic in initial_conditions
                    raw_result =
                        Feather.read(variable_ref)[end, Symbol(PSI.device_name(ic))] # last value of last hour
                    initial_cond = value(PSI.get_value(ic))
                    @test isapprox(raw_result, initial_cond; atol = 1e-2)
                end
            end
        end

        @testset "Verify simulation events" begin
            file = joinpath(
                file_path,
                PSI.get_name(sim),
                output_dir,
                "recorder",
                "simulation.log",
            )
            @test isfile(file)
            events = PSI.list_simulation_events(
                PSI.InitialConditionUpdateEvent,
                joinpath(file_path, "aggregation", output_dir);
                step = 1,
            )
            @test length(events) == 0
            events = PSI.list_simulation_events(
                PSI.InitialConditionUpdateEvent,
                joinpath(file_path, "aggregation", output_dir);
                step = 2,
            )
            @test length(events) == 10
            PSI.show_simulation_events(
                devnull,
                PSI.InitialConditionUpdateEvent,
                joinpath(file_path, "aggregation", output_dir);
                step = 2,
            )
            events = PSI.list_simulation_events(
                PSI.InitialConditionUpdateEvent,
                joinpath(file_path, "aggregation", output_dir);
                step = 1,
                stage = 1,
            )
            @test length(events) == 0
            events = PSI.list_simulation_events(
                PSI.InitialConditionUpdateEvent,
                joinpath(file_path, "aggregation", output_dir);
                step = 2,
                stage = 1,
            )
            @test length(events) == 10
            PSI.show_simulation_events(
                devnull,
                PSI.InitialConditionUpdateEvent,
                joinpath(file_path, "aggregation", output_dir);
                step = 2,
                stage = 1,
            )
        end
    end

    ####################
    stages_definition = Dict(
        "UC" => Stage(
            GenericOpProblem,
            template_hydro_basic_uc,
            c_sys5_hy_uc,
            GLPK_optimizer,
        ),
        "ED" =>
            Stage(GenericOpProblem, template_hydro_ed, c_sys5_hy_ed, ipopt_optimizer),
    )

    sequence = SimulationSequence(
        order = Dict(1 => "UC", 2 => "ED"),
        step_resolution = Hour(1),
        feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon()),
        horizons = Dict("UC" => 24, "ED" => 12),
        intervals = Dict(
            "UC" => (Hour(1), RecedingHorizon()),
            "ED" => (Minute(5), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_from_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterStageChronology(),
    )

    sim = Simulation(
        name = "receding_results",
        steps = 2,
        stages = stages_definition,
        stages_sequence = sequence,
        simulation_folder = file_path,
    )
    build!(sim)
    sim_results = execute!(sim)

    @testset "Test verify time gap for Receding Horizon" begin
        names = ["UC"] # TODO why doesn't this work for ED??
        for name in names
            variable_list = PSI.get_variable_names(sim, name)
            reference_1 = PSI.get_reference(sim_results, name, 1, variable_list[1])[1]
            reference_2 = PSI.get_reference(sim_results, name, 2, variable_list[1])[1]
            time_file_path_1 = joinpath(dirname(reference_1), "time_stamp.feather") #first line, file path
            time_file_path_2 = joinpath(dirname(reference_2), "time_stamp.feather")
            time_1 = convert(Dates.DateTime, Feather.read(time_file_path_1)[1, 1]) # first time
            time_2 = convert(Dates.DateTime, Feather.read(time_file_path_2)[1, 1])
            time_change = time_2 - time_1
            interval = PSI.get_stage_interval(PSI.get_sequence(sim), name)
            @test Dates.Hour(time_change) == Dates.Hour(interval)
        end
    end

    @testset "Test verify parameter feedforward for Receding Horizon" begin
        P_keys = [(PSI.ON, PSY.ThermalStandard)]
        vars_names = [PSI.variable_name(PSI.ON, PSY.ThermalStandard)]
        for (ik, key) in enumerate(P_keys)
            variable_ref = PSI.get_reference(sim_results, "UC", 2, vars_names[ik])[1]
            raw_result = Feather.read(variable_ref)
            ic = PSI.get_parameter_array(PSI.get_parameter_container(
                sim.stages["ED"].internal.psi_container,
                Symbol(key[1]),
                key[2],
            ))
            for name in DataFrames.names(raw_result)
                result = raw_result[1, name] # first time period of results  [time, device]
                initial = value(ic[String(name)]) # [device, time]
                @test isapprox(initial, result, atol = 1.0e-4)
            end
        end
    end

    @testset "Test verify initial condition feedforward for Receding Horizon" begin
        results = load_simulation_results(sim_results, "ED")
        ic_keys = [PSI.ICKey(PSI.DevicePower, PSY.ThermalStandard)]
        vars_names = [PSI.variable_name(PSI.ACTIVE_POWER, PSY.ThermalStandard)]
        ed_horizon = PSI.get_stage_horizon(sim.sequence, "ED")
        no_steps = PSI.get_steps(sim)
        for (ik, key) in enumerate(ic_keys)
            initial_conditions =
                get_initial_conditions(PSI.get_psi_container(sim, "UC"), key)
            vars = results.variable_values[vars_names[ik]] # change to getter function
            for ic in initial_conditions
                output = vars[ed_horizon * (no_steps - 1), Symbol(PSI.device_name(ic))] # change to getter function
                initial_cond = value(PSI.get_value(ic))
                @test isapprox(output, initial_cond, atol = 1.0e-4)
            end
        end
    end
    @testset "Test print methods" begin
        list = [sim, sim_results, sim.sequence, sim_single.sequence, res, sim.stages["UC"]]
        _test_plain_print_methods(list)
        _test_html_print_methods([res])
    end
    @testset "Test print methods of sequence ascii art" begin

        sequence_2 = SimulationSequence(
            order = Dict(1 => "UC", 2 => "ED"),
            step_resolution = Hour(1),
            feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 2)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(1), RecedingHorizon()),
                "ED" => (Minute(5), RecedingHorizon()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_from_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )

        sequence_4 = SimulationSequence(
            order = Dict(1 => "UC", 2 => "ED"),
            step_resolution = Hour(1),
            feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 4)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(1), RecedingHorizon()),
                "ED" => (Minute(5), RecedingHorizon()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_from_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )

        sequence_3 = SimulationSequence(
            order = Dict(1 => "UC", 2 => "ED"),
            step_resolution = Hour(1),
            feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 3)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(1), RecedingHorizon()),
                "ED" => (Minute(5), RecedingHorizon()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_from_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )

        sequence_5 = SimulationSequence(
            order = Dict(1 => "UC", 2 => "ED"),
            step_resolution = Hour(1),
            feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 2)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(1), RecedingHorizon()),
                "ED" => (Minute(5), RecedingHorizon()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => RangeFF(
                    variable_from_stage_ub = PSI.ON,
                    variable_from_stage_lb = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )

        sequence_13 = SimulationSequence(
            order = Dict(1 => "UC", 2 => "ED"),
            step_resolution = Hour(1),
            feedforward_chronologies = Dict(
                ("UC" => "ED") => RecedingHorizon(periods = 13),
            ),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(1), RecedingHorizon()),
                "ED" => (Minute(5), RecedingHorizon()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_from_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        list = [sequence_2, sequence_3, sequence_4, sequence_5, sequence_13]
        _test_plain_print_methods(list)
        stage_1 = FakeStagesStruct(Dict(1 => 1, 2 => 2, 3 => 3, 4 => 4, 5 => 5, 6 => 6)) # testing 5 stages
        stage_3 = FakeStagesStruct(Dict(1 => 1, 2 => 100)) #testing 3 digits
        stage_4 = FakeStagesStruct(Dict(1 => 1, 2 => 1000)) #testing 4 digits
        stage_12 = FakeStagesStruct(Dict(1 => 1, 2 => 12, 3 => 5, 4 => 6))
        list = [stage_1, stage_3, stage_4, stage_12]
        _test_plain_print_methods(list)

    end

    ####################
    @testset "negative test checking total sums" begin
        stage_names = keys(sim.stages)
        for name in stage_names
            files = collect(readdir(sim_results.results_folder))
            for f in files
                rm("$(sim_results.results_folder)/$f")
            end
            variable_list = PSI.get_variable_names(sim, name)
            res = load_simulation_results(sim_results, name)
            write_results(res)
            _file_path = joinpath(sim_results.results_folder, "$(variable_list[1]).feather")
            rm(_file_path)
            fake_df = DataFrames.DataFrame(:A => Array(1:10))
            Feather.write(_file_path, fake_df)
            @test_logs(
                (:error, r"hash mismatch"),
                match_mode = :any,
                @test_throws(
                    IS.HashMismatchError,
                    check_file_integrity(dirname(_file_path))
                )
            )
        end
        for name in stage_names
            variable_list = PSI.get_variable_names(sim, name)
            check_file_path = PSI.get_reference(sim_results, name, 1, variable_list[1])[1]
            rm(check_file_path)
            time_length = sim_results.chronologies["stage-$name"]
            fake_df = DataFrames.DataFrame(:A => Array(1:time_length))
            Feather.write(check_file_path, fake_df)
            @test_logs(
                (:error, r"hash mismatch"),
                match_mode = :any,
                @test_throws(
                    IS.HashMismatchError,
                    check_file_integrity(dirname(check_file_path))
                )
            )
        end
    end

    @testset "Simulation with Cache" begin
        stages_definition = Dict(
            "UC" => Stage(
                GenericOpProblem,
                template_hydro_st_standard_uc,
                #template_uc,
                c_sys5_hy_uc,
                #c_sys5_uc,
                GLPK_optimizer,
            ),
            "ED" => Stage(
                GenericOpProblem,
                template_hydro_st_ed,
                #template_ed,
                c_sys5_hy_ed,
                #c_sys5_ed,
                GLPK_optimizer,
            ),
        )

        sequence_cache = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_from_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
                ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                    variable_from_stage = PSI.ACTIVE_POWER,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            cache = Dict(
                ("UC",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON),
                ("UC", "ED") => StoredEnergy(PSY.HydroEnergyReservoir, PSI.ENERGY),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim_cache = Simulation(
            name = "cache",
            steps = 2,
            stages = stages_definition,
            stages_sequence = sequence_cache,
            simulation_folder = file_path,
        )
        build!(sim_cache)
        sim_cache_results = execute!(sim_cache)
        var_names =
            axes(PSI.get_stage(sim_cache, "UC").internal.psi_container.variables[:On__ThermalStandard])[1]
        for name in var_names
            var =
                PSI.get_stage(sim_cache, "UC").internal.psi_container.variables[:On__ThermalStandard][
                    name,
                    24,
                ]
            cache = PSI.get_cache(
                sim_cache,
                PSI.CacheKey(TimeStatusChange, PSY.ThermalStandard),
            ).value[name]
            @test JuMP.value(var) == cache[:status]
        end

        @testset "Test verify initial condition update using StoredEnergy cache" begin
            ic_keys = [PSI.ICKey(PSI.EnergyLevel, PSY.HydroEnergyReservoir)]
            vars_names = [PSI.variable_name(PSI.ENERGY, PSY.HydroEnergyReservoir)]
            for (ik, key) in enumerate(ic_keys)
                variable_ref =
                    PSI.get_reference(sim_cache_results, "ED", 1, vars_names[ik])[end]
                initial_conditions =
                    get_initial_conditions(PSI.get_psi_container(sim_cache, "UC"), key)
                for ic in initial_conditions
                    raw_result =
                        Feather.read(variable_ref)[end, Symbol(PSI.device_name(ic))] # last value of last hour
                    initial_cond = value(PSI.get_value(ic))
                    @test isapprox(raw_result, initial_cond)
                end
            end
        end
    end

    @testset "" begin
        single_stage_definition = Dict(
            "ED" => Stage(
                GenericOpProblem,
                template_hydro_st_ed,
                c_sys5_hy_ed,
                GLPK_optimizer,
            ),
        )

        single_sequence = SimulationSequence(
            step_resolution = Hour(1),
            order = Dict(1 => "ED"),
            horizons = Dict("ED" => 12),
            intervals = Dict("ED" => (Hour(1), Consecutive())),
            cache = Dict(("ED",) => StoredEnergy(PSY.HydroEnergyReservoir, PSI.ENERGY)),
            ini_cond_chronology = IntraStageChronology(),
        )

        sim_single = Simulation(
            name = "cache_st",
            steps = 2,
            stages = single_stage_definition,
            stages_sequence = single_sequence,
            simulation_folder = file_path,
        )
        build!(sim_single)
        sim_cache_results = execute!(sim_single)

        @testset "Test verify initial condition update using StoredEnergy cache" begin
            ic_keys = [PSI.ICKey(PSI.EnergyLevel, PSY.HydroEnergyReservoir)]
            vars_names = [PSI.variable_name(PSI.ENERGY, PSY.HydroEnergyReservoir)]
            for (ik, key) in enumerate(ic_keys)
                variable_ref =
                    PSI.get_reference(sim_cache_results, "ED", 1, vars_names[ik])[1]
                initial_conditions =
                    get_initial_conditions(PSI.get_psi_container(sim_single, "ED"), key)
                for ic in initial_conditions
                    raw_result =
                        Feather.read(variable_ref)[end, Symbol(PSI.device_name(ic))] # last value of last hour
                    initial_cond = value(PSI.get_value(ic))
                    @test isapprox(raw_result, initial_cond)
                end
            end
        end
    end
end

@testset "Test load simulation" begin
    path = (joinpath(pwd(), "test_reading_results"))
    !isdir(path) && mkdir(path)

    try
        test_load_simulation(path)
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
using DEBTrait
using CSV, DataFrames, Statistics
using DifferentialEquations, Plots

include("/Users/glmarschmann/.julia/dev/DEBTrait/parameterization/positive_negative.jl")

u0 = zeros(p.dim)
u0[1+p.n_polymers:p.n_polymers+p.n_monomers] = 25.0*ones(p.n_monomers)
u0[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes] = p.Setup.u0_reserves
u0[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes] = p.Setup.u0_structures

p.PMetabolismC.α = [0.01, 0.3]
p.PMetabolismC.y_EX = [1.4, 2.8]

tspan = (0.0,3.3)
prob  = ODEProblem(rhs_BatchC!,u0,tspan,p)
sol   = solve(prob, alg_hints=[:stiff])

plot(sol, vars=1+p.n_polymers:p.n_polymers+p.n_monomers)
plot(sol, vars=1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes)
plot(sol, vars=1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes)
plot(sol, vars=1+p.n_polymers+p.n_monomers+2*p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes)
plot(sol, vars=1+p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes+p.n_microbes)

BGE_pos = [DEBTrait.calc_BGE(zeros(p.dim), sol.u[i], p, 0.0)[1] for i in 1:size(sol.t,1)]
BGE_neg = [DEBTrait.calc_BGE(zeros(p.dim), sol.u[i], p, 0.0)[2] for i in 1:size(sol.t,1)]
BGE_median = [median(BGE_pos), median(BGE_neg)]
plot(sol.t, BGE_pos)
plot!(sol.t, BGE_neg)

r_pos = [DEBTrait.calc_growth_rate(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
r_neg = [DEBTrait.calc_growth_rate(sol.u[i], p)[2] for i in 1:size(sol.t,1)]
r_median = [median(r_pos), median(r_neg)]
plot(sol.t, r_pos)
plot!(sol.t, r_neg)

scatter(r_median, BGE_median)

x_pos = [DEBTrait.calc_growth_production(sol.u[i], p)[1][1] for i in 1:size(sol.t,1)]
x_neg = [DEBTrait.calc_growth_production(sol.u[i], p)[1][2] for i in 1:size(sol.t,1)]
plot(sol.t, x_pos)
plot!(sol.t, x_neg)

ext: jl


text:
struct CombinedConfusionMatrix
    matrices
    classlabels

    function CombinedConfusionMatrix(matrices;
                                                    classlabels=nothing)
        if classlabels == nothing
            _, n = size(matrices[1])
            classlabels = ["class $x" for x in 1:n]
        end

        new(matrices, classlabels)
    end

end

function checkedscore(a,b)
    b==0? 0:a/b
end

function score(c)
    m,n = size(c)
    sum([c[i,i] for i in 1:m])
end

function expandmatrix(matrix::Matrix)

    m, n = size(matrix)

    result = Matrix{Number}(m+2, n+2)

    result .= 0

    for row in 1 : m
        for col in 1: n
            result[row,col] = matrix[row,col]
        end

        result[row, n+1] = sum(matrix[row,:])
        rowscore = checkedscore(matrix[row,row], sum(matrix[row,:]))
        result[row, n+2] = rowscore
    end

    for col in 1: n
        result[m+1, col] = sum(matrix[:,col])
        colscore = checkedscore( matrix[col,col], sum(matrix[:,col]))
        result[m+2, col] = colscore
    end

    result[m+1, n+1] = sum(matrix)

    totalscore =   score(matrix) / sum(matrix)
    result[m+2, n+2] = totalscore

    result

end

function tohtml(confusionmatrix::CombinedConfusionMatrix)

    cms = confusionmatrix.matrices
    labels = confusionmatrix.classlabels

    cm = cms[1]
    m,n = size(cms[1])

    cms2 = expandmatrix.(cms)
    d = Float64.(cat(3, cms2...))
    sems = mapslices(StatsBase.sem, d, 3)
    means = mapslices(StatsBase.mean, d, 3)

    stream = IOBuffer()
    write(stream,"""
""")
    write(stream,"
")
    write(stream,"""	""")
    write(stream,"""Predicted class x&#770	""")
    write(stream,"
")
    write(stream,"
")
    for label in labels
        write(stream,"$(label)	")
    end
    write(stream,"Row total	")
    write(stream,"Accuracy  P(x&#770|x)	")
    write(stream,"
")
    write(stream,"
")
    write(stream,"""""")
    write(stream,"""
Actual

class x
""")
    write(stream,"	")
    write(stream,"$(labels[1])	")
    for j in 1:n
        write(stream,"""$(@sprintf("%.2f",means[1,j])) ± $(@sprintf("%.2f",sems[1,j]))	""")
    end
    write(stream,"""$(@sprintf("%.2f",means[1,n+1])) ± $(@sprintf("%.2f",sems[1,n+1]))	""") # Row total
    write(stream,"""$(@sprintf("%.2f",means[1,n+2])) ± $(@sprintf("%.2f",sems[1,n+2]))	""") # Accuracy
    write(stream,"
")
    for i in 2:length(labels)
        write(stream,"
")
        write(stream,"$(labels[i])	")
        for j in 1:n
            write(stream,"""$(@sprintf("%.2f",means[i,j])) ± $(@sprintf("%.2f",sems[i,j]))	""")
        end
        write(stream,"""$(@sprintf("%.2f",means[i,n+1])) ± $(@sprintf("%.2f",sems[i,n+1]))	""") # Row total
        write(stream,"""$(@sprintf("%.2f",means[i,n+2])) ± $(@sprintf("%.2f",sems[i,n+2]))	""") # Accuracy
        write(stream,"
")
    end
    write(stream,"
")
    write(stream,"""	""")
    write(stream,"Column total	")
    for j in 1:length(labels)
        write(stream,"""$(@sprintf("%.2f",means[m+1,j])) ± $(@sprintf("%.2f",sems[m+1,j]))	""")
    end
    write(stream,"""$(@sprintf("%.2f",means[m+1,n+1])) ± $(@sprintf("%.2f",sems[m+1,n+1]))	""")
    write(stream,"
")
    write(stream,"
")
    write(stream,"""	""")
    write(stream,"Reliability P(x|x&#770)	")
    for j in 1:length(labels)
        write(stream,"""$(@sprintf("%.2f",means[m+2,j])) ± $(@sprintf("%.2f",sems[m+2,j]))	""") # Reliability
    end
    write(stream,"""	""")
    write(stream,"""$(@sprintf("%.2f",means[m+2,n+2])) ± $(@sprintf("%.2f",sems[m+2,n+2]))	""")
    write(stream,"
")
    write(stream,"
")
    String(stream)
end

Base.show(io::IO, ::MIME"text/html", c::CombinedConfusionMatrix) =
    print(io, tohtml(c))


text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

using SpecialFunctions
using ChainRulesCore
using ChainRulesTestUtils
using Random
using Test
using Base.MathConstants: γ

using SpecialFunctions: AmosException, f64

# useful test functions for relative error, which differ from isapprox
# relerr separately looks at the real and imaginary parts if one of the arguments is complex
relerr(z, x) = z == x ? zero(z) : abs(z - x) / abs(x)
relerr(z::T, x::T) where {T <: Complex} = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
relerr(z::Complex, x) = relerr(z, Complex(x))
relerr(z, x::Complex) = relerr(Complex(z), x)
checktol(err::Float16) = err ≤ 5e-2
checktol(err::Float32) = err ≤ 1e-6
checktol(err::Float64) = err ≤ 1e-13
≅(a,b) = checktol(relerr(a,b))

tests = [
    "bessel",
    "beta_inc",
    "betanc",
    "ellip",
    "erf",
    "expint",
    "gamma_inc",
    "gamma",
    "sincosint",
    "other_tests",
    "chainrules"
]

const testdir = dirname(@__FILE__)

for t in tests
    tp = joinpath(testdir, "$(t).jl")
    @testset "$(t)" begin
      include(tp)
    end
end


text:
using Documenter, Lindenmayer

makedocs(
    sitename = "Lindenmayer",
    modules  = [Lindenmayer],
    format   = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages    = Any[
        "Lindenmayer.jl" => "index.md",
        "Examples"       => "examples.md",
        ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo   = "github.com/cormullion/Lindenmayer.jl.git",
    target = "build"
)

max_stars_repo_name: cormullion/Lindenmayer

max_issues_repo_name: cormullion/Lindenmayer

max_forks_repo_name: cormullion/Lindenmayer


text:
# Autogenerated wrapper script for libmicrohttpd_jll for i686-w64-mingw32
export libmicrohttpd

using GnuTLS_jll
JLLWrappers.@generate_wrapper_header("libmicrohttpd")
JLLWrappers.@declare_library_product(libmicrohttpd, "libmicrohttpd-12.dll")
function __init__()
    JLLWrappers.@generate_init_header(GnuTLS_jll)
    JLLWrappers.@init_library_product(
        libmicrohttpd,
        "bin\\libmicrohttpd-12.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m5(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m5!(A_copy)
end

@inline function exp_ps_m5!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=6
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 A3 A4 P2 C1 P1 C0 P0
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A3 with operation: mult
    mul!(memslots3,memslots1,memslots2)
    # Computing A4 with operation: mult
    mul!(memslots4,memslots1,memslots3)
    # Computing P2 = x*I+x*A+x*A2+x*A3+x*A4
    coeff1=2.48015873015873e-5
    coeff2=2.7557319223985893e-6
    coeff3=2.755731922398589e-7
    coeff4=2.505210838544172e-8
    coeff5=2.08767569878681e-9
    memslots5 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4
    mul!(memslots5,true,I*coeff1,true,true)
    # Computing C1 with operation: mult
    mul!(memslots6,memslots5,memslots4)
    # Deallocating P2 in slot 5
    # Computing P1 = x*I+x*A+x*A2+x*A3+x*C1
    coeff1=0.041666666666666664
    coeff2=0.008333333333333333
    coeff3=0.001388888888888889
    coeff4=0.0001984126984126984
    coeff5=1.0
    # Smart lincomb recycle C1
    memslots6 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots6
    mul!(memslots6,true,I*coeff1,true,true)
    # Computing C0 with operation: mult
    mul!(memslots5,memslots6,memslots4)
    # Deallocating P1 in slot 6
    # Deallocating A4 in slot 4
    # Computing P0 = x*I+x*A+x*A2+x*A3+x*C0
    coeff1=1.0
    coeff2=1.0
    coeff3=0.5
    coeff4=0.16666666666666666
    coeff5=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots5
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating A2 in slot 2
    # Deallocating A3 in slot 3
    # Deallocating C0 in slot 5
    return memslots1 # Returning P0
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using PyQHA
using Documenter

DocMeta.setdocmeta!(PyQHA, :DocTestSetup, :(using PyQHA); recursive=true)

makedocs(;
    modules=[PyQHA],
    authors="Qi Zhang <singularitti@outlook.com>",
    repo="https://github.com/MineralsCloud/PyQHA.jl/blob/{commit}{path}#{line}",
    sitename="PyQHA.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MineralsCloud.github.io/PyQHA.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MineralsCloud/PyQHA.jl",
)


text:
# Insertion sort
function insertionsort!(arr, low, high)
    @inbounds for i = (low+1):high
        j = i
        x = arr[i]
        while j > low && arr[j] < arr[j-1]
            exch!(arr, j, j-1)
            j -= 1
        end
    end
end

insertsort!(arr) = insertionsort!(arr, 1, length(arr))
const SMALL_THRESHOLD = 20

# Merge sort
function mergesort!(arr)
    mergesort!(arr, 1, length(arr))
end

function mergesort!(arr, low, high)
    if low < high
        middle = Int(floor((low + high) / 2))
        mergesort!(arr, low, middle)
        mergesort!(arr, middle+1, high)
        merge!(arr, low, middle, high)
    end
end

function merge!(arr, low, middle, high)
    leftarr = copy(arr[low:middle])
    rightarr = copy(arr[(middle+1):high])

    i = low
    # compare the first elements between two arrays
    while (!(isempty(leftarr) || isempty(rightarr)))
        if leftarr[1] <= rightarr[1]
            arr[i] = popfirst!(leftarr)
        else
            arr[i] = popfirst!(rightarr)
        end
        i += 1
    end

    # add remaining elements
    while !isempty(leftarr)
        arr[i] = popfirst!(leftarr)
        i += 1
    end
    while !isempty(rightarr)
        arr[i] = popfirst!(rightarr)
        i += 1
    end
end

"""
    quicksort!(arr)

Quicksort implemntations based on Algorithm Design Manual and Julia's
Official implementation.

This implementation is about 10% slow than Julia's `sort!`
(Quicksort + Insertionsort). The pivot selecting mechanism follows Julia's
implementaiton that select three points and pick up the middle one, instead of
shuffling the entire list. One bottleneck is in Julia they have implemented
a faster float comparison: split the list into negative and positve, and then
sort with `slt_int`.

See [my gist](https://gist.github.com/Jetafull/977eaf8a501891b2527a950d777efb5a)
for more details.
"""
function quicksort!(arr::AbstractVector)
    quicksort!(arr, 1, length(arr))
end

# standard implementation:
# pivot is the final position, recursively sort on both sides of the pivot.
function quicksort!(arr, lo, hi)
    lenarr = hi - lo + 1
    if lenarr <= SMALL_THRESHOLD
        insertionsort!(arr, lo, hi)
    else
        if lo < hi
            firsthigh = partition!(arr, lo, hi)
            quicksort!(arr, lo, firsthigh-1)
            quicksort!(arr, firsthigh+1, hi)
        end
    end
end

lt(x, y) = x < y

@inline function exch!(arr, i, j)
    @inbounds begin
        arr[i], arr[j] = arr[j], arr[i]
    end
end

# Baseline version
# Algorithm Design Manual: 4.6
# Need to shuffle the entire list before sorting (uncomment shuffle!(arr))
# function partition!(arr, low, high)
#     pivot = high
#     firsthigh = low
#     i = low
#     @inbounds while i < high
#         if arr[i] < arr[pivot]
#             exch!(arr, i, firsthigh)
#             firsthigh += 1
#         end
#         i += 1
#     end
#     exch!(arr, pivot, firsthigh)
#     return firsthigh
# end

# Julia offical implementation:
# https://github.com/JuliaLang/julia/blob/master/base/sort.jl
function partition!(v, lo, hi)
    @inbounds begin
        pivot = v[selectpivot!(v, lo, hi)]
        i, j = lo, hi
        while true
            i += 1; j -= 1
            while lt(v[i], pivot); i += 1; end;
            while lt(pivot, v[j]); j -= 1; end;
            i >= j && break
            v[i], v[j] = v[j], v[i]
        end
        v[j], v[lo] = pivot, v[j]
        return j
    end
end

midpoint(lo, hi) = lo + ((hi - lo) >>> 0x01)
# Given 3 locations in an array (lo, mi, and hi), sort v[lo], v[mi], v[hi]) and
# choose the middle value as a pivot
#
# Upon return, the pivot is in v[lo], and v[hi] is guaranteed to be
# greater than the pivot
@inline function selectpivot!(v, lo, hi)
    @inbounds begin
        mi = midpoint(lo, hi)
        # sort v[mi] <= v[lo] <= v[hi] such that the pivot is immediately in place
        if lt(v[lo], v[mi])
            v[mi], v[lo] = v[lo], v[mi]
        end

        if lt(v[hi], v[lo]) # v[lo] > v[hi] and v[lo] > v[mi]
            if lt(v[hi], v[mi])
                v[hi], v[lo], v[mi] = v[lo], v[mi], v[hi]
            else
                v[hi], v[lo] = v[lo], v[hi]
            end
        end

        return lo
    end
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
import Base.Dict

export Atoms

"""
`JData`: datatype for storing any data

some data which needs to be updated if the configuration (positions only!) has
changed too much.
"""
mutable struct JData{T}
   max_change::T     # how much X may change before recomputing
   accum_change::T   # how much has it changed already
   data::Any
end

struct LinearConstraint{T}
   C::Matrix{T}
   b::Vector{T}
   desc::String   # description of the constraint
end

mutable struct DofManager{T}
   variablecell::Bool
   xfree::Vector{Int}
   lincons::Vector{LinearConstraint{T}}
   # ----
   X0::Vector{JVec{T}}
   F0::JMat{T}
end

"""
`Atoms{T <: AbstractFloat} <: AbstractAtoms`

The main type storing information about atomic configurations. This is
not normally constructed directly by calling `Atoms`, but using one of the
utility functions such as `bulk`.

## Constructors

TODO

## Convenience functions

* `bulk`

# Getters and Setters

* `positions`, `set_positions!`
* `momenta`, `set_momenta!`
* `masses`, `set_masses!`
* `numbers`, `set_numbers!`
* `cell`, `set_cell!`
* `pbc`, `set_pbc!`
* `calculator`, `set_calculator!`
* `get_data`, `set_data!`
"""
mutable struct Atoms{T <: AbstractFloat} <: AbstractAtoms{T}
   X::Vector{JVec{T}}              # positions
   P::Vector{JVec{T}}              # momenta (or velocities?)
   M::Vector{T}                    # masses
   Z::Vector{AtomicNumber}         # atomic numbers
   cell::JMat{T}                   # cell
   pbc::JVec{Bool}                 # boundary condition
   calc::Union{Nothing, AbstractCalculator}
   dofmgr::DofManager
   data::Dict{Any,JData{T}}
end

Atoms(X::Vector{JVec{T}}, P::Vector{JVec{T}}, M::Vector{T},
      Z::Vector{AtomicNumber}, cell::JMat{T}, pbc::JVec{Bool},
      calc::Union{Nothing, AbstractCalculator} = nothing) where {T} =
   Atoms(X, P, M, Z, cell, pbc, calc,
         DofManager(length(X), T), Dict{Any,JData{T}}())

fltype(::Atoms{T}) where {T} = T

# derived properties
length(at::Atoms) = length(at.X)

getindex(at::Atoms, i::Integer) = at.X[i]

function setindex!(at::Atoms{T}, val, i::Integer) where {T}
   u = norm(at.X[i] - val)
   at.X[i] = convert(JVec{T}, val)
   update_data!(at, u)
   return val
end

# ------------------------ access to struct fields ----------------------

# ----- getters

positions(at::Atoms) = copy(at.X)
momenta(at::Atoms) = copy(at.P)
masses(at::Atoms) = copy(at.M)
numbers(at::Atoms) = copy(at.Z)
atomic_numbers(at::Atoms) = copy(at.Z)
cell(at::Atoms) = at.cell
pbc(at::Atoms) = at.pbc
calculator(at::Atoms) = at.calc

chemical_symbols(at::Atoms) = Chemistry.chemical_symbol.(at.Z)

# ----- setters

function set_positions!(at::Atoms{T}, X::AbstractVector{JVec{T}})  where T
   update_data!(at, dist(at, X))
   at.X .= X
   return at
end
function set_momenta!(at::Atoms{T}, P::AbstractVector{JVec{T}})  where T
   at.P .= P
   return at
end
function set_masses!(at::Atoms{T}, M::AbstractVector{T})  where T
   update_data!(at, Inf)
   at.M .= M
   return at
end
function set_numbers!(at::Atoms, Z::AbstractVector{<:Integer})
   update_data!(at, Inf)
   at.Z .= AtomicNumber.(Z)
   return at
end
function set_cell!(at::Atoms{T}, C::AbstractMatrix) where T
   update_data!(at, Inf)
   at.cell = JMat{T}(C)
   return at
end
function set_pbc!(at::Atoms, p::Union{AbstractVector, Tuple})
   q = JVec{Bool}(p...)
   if at.pbc != q
      update_data!(at, Inf)
      at.pbc = q
   end
   return at
end
function set_calculator!(at::Atoms, calc::Union{AbstractCalculator, Nothing})
   at.calc = calc
   return at
end

# --------------- equality tests -------------------

import Base.==
==(at1::Atoms{T}, at2::Atoms{T}) where T = (
   isapprox(at1, at2, tol = 0.0)  &&
          (at1.data == at2.data)  &&
          (at1.calc == at2.calc)  &&
          (at1.dofmgr == at2.dofmgr) )

import Base.isapprox
function isapprox(at1::Atoms{T}, at2::Atoms{T}; tol = sqrt(eps(T))) where {T}
   if length(at1) != length(at2)
      return false
   end
   if tol > 0
      ndigits = floor(Int, abs(log10(tol)))
      X1 = [round.(x; digits=ndigits) for x in at1.X]
      X2 = [round.(x; digits=ndigits) for x in at2.X]
   else
      X1, X2 = at1.X, at2.X
   end
   p1 = sortperm(X1)
   p2 = sortperm(X2)
   return (maxdist(at1.X[p1], at2.X[p2]) <= tol) &&
          (maxdist(at1.P[p1], at2.P[p2]) <= tol)  &&
          (norm(at1.M[p1] - at2.M[p2], Inf) <= tol) &&
          (at1.Z[p1] == at2.Z[p2]) &&
          (norm(at1.cell - at2.cell, Inf) <= tol) &&
          (at1.pbc == at2.pbc)
end

# --------------- some aliases for easier / direct access --------------
#                 and some access not covered above

function neighbourlist(at::Atoms{T}, rcut::AbstractFloat;
            recompute=false, key="nlist:default", storelist=true,
            int_type = Int, kwargs...) where {T}
   # we are allowed to use a stored list
   if !recompute
      if has_data(at, key)
         nlist = get_data(at, key)::PairList{T, int_type}
         if cutoff(nlist) >= rcut
            return nlist
         end
      end
   end
   # we are either forces to recompute, or no nlist exists or is not
   # suitable for the request.
   nlist = PairList(positions(at), rcut, cell(at), pbc(at);
                    int_type=int_type, kwargs...)::PairList{T, int_type}
   if storelist
      set_data!(at, key, nlist, 0.0)
   end
   return nlist
end

# function neighbourlist(at::Atoms{T}, rcut::AbstractFloat;
#             recompute=false, key="nlist:default", kwargs...) where {T}
#    return PairList(positions(at), rcut, cell(at), pbc(at); kwargs...)
# end

neighbourlist(at::Atoms) = neighbourlist(at, cutoff(at))

# --------------- data Dict handling ------------------

"""
`has_data(at::Atoms, key)`: checks whether `at` is storing a datum with key `key`; see also `set_data!`
"""
has_data(at::Atoms, key) = haskey(at.data, key)

"""
`get_data(at::Atoms, key)`: retrieves a datum;
see also `set_data!`. Unlike `get_positions`, etc, this returns the original
datum, not a copy (unless it is immutable).
"""
get_data(at::Atoms, key) = (at.data[key]).data

"""
`set_data!(a::Atoms, key, value, max_change=Inf)`:

Stores `value` under an arbitrary key `key` inside `at`. The entry is paired
with an updatemeasure `chg`. Whenever atom positions or the cell change by a
distance `d`, the counter is incremented `chg += d`. When `chg > max_change`,
the entry is deleted from the dictionary.

The two most common uses are with `max_change = Inf` (default) - in this case
the entry will never be deleted; and `max_change = 0.0` - in this case the entry
will be deleted whenever atom positions or the cell change.

The change counter is used, e.g., for neighbourlists with buffer, and for
preconditioners so that the list need not be recomputed each time the
configuration changes.
"""
function set_data!(at::Atoms, key, value, max_change=Inf)
   at.data[key] = JData(max_change, 0.0, value)
   return at
end

"""
`function set_transient!(a::Atoms, key, value, max_change=0.0)`: this is an
alias for `set_data!` but with default `max_change = 0.0` instead of `max_change =
Inf`.
"""
set_transient!(at::Atoms, key, value, max_change=0.0) =
      set_data!(at, key, value, max_change)

"""
`update_data!(at::Atoms, r::Real)`:
increment the change counter in all stored data and delete the
expired ones.
"""
function update_data!(at::Atoms, r::Real)
   for (key, t) in at.data
      t.accum_change += r
      if t.accum_change >= t.max_change
         delete!(at.data, key)
      end
   end
   return at
end

# ------------------------ FIO  ----------------------

JuLIP.FIO.write_dict(at::Atoms) =
   Dict( "__id__" => "JuLIP_Atoms",
         "X"      => at.X,
         "P"      => at.P,
         "M"      => at.M,
         "Z"      => Int.(at.Z),
         "cell"   => at.cell,
         "pbc"    => at.pbc,
         "calc"   => nothing,
         "data"   => nothing )

Atoms(D::Dict) = Atoms(D["X"], D["P"], D["M"], AtomicNumber.(D["Z"]),
                       D["cell"], D["pbc"])
          # calc = D["calc"],
          # cons = D["cons"],
          # data = D["data"] )

JuLIP.FIO.read_dict(::Val{:JuLIP_Atoms}, D::Dict) = Atoms(D)

import ExtXYZ

# conversion rules from Atoms.data entries to ExtXYZ standard types
_write_convert(value) = nothing # default rule: do not write
_write_convert(value::AbstractVector) = value
_write_convert(value::AbstractVector{JVec{T}}) where T = value |> mat |> Matrix
_write_convert(value::JVec{T}) where T = value |> Array
_write_convert(value::JMat{T}) where T = value |> Matrix
_write_convert(value::Union{Number,String}) = value

# (value, natoms) -> destination ∈ (:info, :arrays, missing)
_dest(::Any, natoms::Int) = missing

# this is a little fragile, since it assumes all vectors of length `natoms` contain per-atom data
# this is problematic for e.g. `natoms=3`.
_dest(value::AbstractVector, natoms::Int) = length(value) == natoms ? :arrays : :info
_dest(::JVec{T}, natoms::Int) where T = :info
_dest(::JMat{T}, natoms::Int) where T = :info
_dest(::Union{Number,String}, natoms::Int) = :info

function _atoms_to_extxyz_dict(atoms::Atoms{T}) where {T}
   dict = write_dict(atoms)
   natoms = length(atoms)
   dict["N_atoms"] = natoms
   dict["cell"] = _write_convert(dict["cell"])
   dict["pbc"] = _write_convert(dict["pbc"])

   # try to figure out whether data entries are per-atom or per-config
   info = Dict{String}{Any}()
   arrays = Dict{String}{Any}()
   for key in keys(atoms.data)
      value = get_data(atoms, key)
      dest, value = _dest(value, natoms), _write_convert(value)
      if value === nothing || ismissing(dest)
         continue
      end
      key = string(key) # symbol -> string, if needed
      if dest == :arrays
         push!(arrays, key => value)
      else
         push!(info, key => value)
      end
   end

   # move Atoms fields into `arrays` sub-dict
   # first X and Z, which should always be initialised
   arrays["pos"] = _write_convert(pop!(dict, "X"))
   arrays["species"] = _write_convert(string.(chemical_symbol.(atoms.Z)))
   arrays["Z"] = _write_convert(pop!(dict, "Z"))

   # now M and P, which might not be set
   for (source, dest) in zip(["M", "P"], ["masses", "momenta"])
      if source ∈ keys(dict) && dict[source] !== nothing && length(dict[source]) == natoms
         arrays[dest] = _write_convert(pop!(dict, source))
      end
   end

   # tidy up top-level dict entries
   dict["info"] = info
   dict["arrays"] = arrays
   delete!(dict, "data") # remove `nothing` entry
   delete!(dict, "calc") # remove `nothing` entry
   delete!(dict, "__id__")

   return dict
end

# define conversions to be applied when reading from ExtXYZ files to JuLIP Atoms
_read_convert(value, ::Int) = value
_read_convert(value::AbstractMatrix, natoms::Int) = size(value, 2) == natoms ? vecs(value) : value

function _extxyz_dict_to_atoms(dict)
   dict["__id__"] = "JuLIP_Atoms"

   natoms = Int(pop!(dict, "N_atoms"))
   info = pop!(dict, "info")
   arrays = pop!(dict, "arrays")

   "pos" in keys(arrays) || error("arrays dictionary missing 'pos' entry containing positions")
   dict["X"] = _read_convert(pop!(arrays, "pos"), natoms)
   @assert length(dict["X"]) == natoms

   # atomic numbers and symbols
   dict["Z"] = "Z" in keys(arrays) ? pop!(arrays, "Z") : nothing
   if "species" in keys(arrays)
      species = _read_convert(pop!(arrays, "species"), natoms)
      Zsp = atomic_number.([Symbol(sp) for sp in species])
      if dict["Z"] !== nothing
         all(dict["Z"] .== Zsp) || error("inconsistent 'Z' and 'species' properties")
      else
         dict["Z"] = Zsp
      end
   end
   dict["Z"] === nothing && error("atomic numbers not defined - either 'Z' or 'species' must be present")

   # mass - lookup from atomic number if not present
   if "masses" in keys(arrays)
      dict["M"] = _read_convert(pop!(arrays, "masses"), natoms)
   else
      dict["M"] = [atomic_mass(z) for z in AtomicNumber.(dict["Z"])]
   end

   # momenta / velocities
   if "momenta" in keys(arrays)
      dict["P"] = _read_convert(pop!(arrays, "momenta"), natoms)
   else
      dict["P"] = _read_convert(zeros((3, nat)), natoms)
   end

   atoms = read_dict(dict)

   # everything else goes in data
   for (label, D) in zip((:info, :arrays), (info, arrays))
      for (key, value) in D
            value = _read_convert(value, natoms)
            set_data!(atoms, key, value)
      end
   end
   return atoms
end

"""
   read_extxyz(file, atoms[, range])

Read from an extended XYZ file using ExtXYZ.jl.

`file` can be a string filename, open IO stream or `FILE*` pointer.
`range` can be absent, an integer frame number, integer range, or array of integers.
Default behaviour is to read all frames in file.

Returns a vector of `Atoms` structs
"""
function read_extxyz(file, args...; kwargs...)
   results = Atoms[]
   for dict in ExtXYZ.read_frames(file, args...; kwargs...)
      push!(results, _extxyz_dict_to_atoms(dict))
   end
   return results
end

"""
   write_extxyz(file, atoms)

Write atoms to file in extended XYZ format using `ExtXYZ.jl`.

`file` can be a string filename, open IO stream or `FILE*` pointer.
`atoms` can be either a single `JuLIP.Atoms` or a vector of `JuLIP.Atoms`.
"""
write_extxyz(file, atoms::Atoms{T}) where T = ExtXYZ.write_frame(file, _atoms_to_extxyz_dict(atoms))

write_extxyz(file, atoms::Vector{Atoms{T}}) where T = ExtXYZ.write_frames(file, _atoms_to_extxyz_dict.(atoms))

export read_extxyz, write_extxyz


text:
using .CommManager

# Don't send previous lines to the completions function,
# due to issue #380.  Find the start of the first line
# (if any) where the expression is parseable.  Replace
# with find_parsestart(c,p) = start(c) once julia#9467 is merged.
parseok(s) = !Meta.isexpr(Meta.parse(s, raise=false), :error)
function find_parsestart(code, cursorpos)
    s = firstindex(code)
    while s < cursorpos
        parseok(code[s:cursorpos]) && return s
        s = nextind(code, s)
        while s < cursorpos && code[s] ∉ ('\n','\r')
            s = nextind(code, s)
        end
    end
    return firstindex(code) # failed to find parseable lines
end

# As described in jupyter/jupyter_client#259, Jupyter's cursor
# positions are actually measured in UTF-16 code units, not
# Unicode characters.  Hence, these functions are similar to
# Base.chr2ind and Base.ind2chr but count non-BMP characters
# as 2 code units.
function utf16_to_ind(str, ic)
    i = 0
    e = lastindex(str)
    while ic > 0 && i < e
        i = nextind(str, i)
        ic -= UInt32(str[i]) < 0x10000 ? 1 : 2
    end
    return i
end
function ind_to_utf16(str, i)
    ic = 0
    i = min(i, lastindex(str))
    while i > 0
        ic += UInt32(str[i]) < 0x10000 ? 1 : 2
        i = prevind(str, i)
    end
    return ic
end

# protocol change in Jupyter 5.2 (jupyter/jupyter_client#262)
chr2ind(m::Msg, str::String, ic::Integer) = ic == 0 ? 0 :
    VersionNumber(m.header["version"]) ≥ v"5.2" ? nextind(str, 0, ic) : utf16_to_ind(str, ic)
ind2chr(m::Msg, str::String, i::Integer) = i == 0 ? 0 :
    VersionNumber(m.header["version"]) ≥ v"5.2" ? length(str, 1, i) : ind_to_utf16(str, i)
#Compact display of types for Jupyterlab completion

import REPL: REPLCompletions
import REPL.REPLCompletions: sorted_keywords, emoji_symbols, latex_symbols

complete_type(::Type{<:Function}) = "function"
complete_type(::Type{<:Type}) = "type"
complete_type(::Type{<:Tuple}) = "tuple"

function complete_type(T::DataType)
    s = string(T)
    (textwidth(s) ≤ 20 || isempty(T.parameters)) && return s
    buf = IOBuffer()
    print(buf, T.name)
    position(buf) > 19 && return String(take!(buf))
    print(buf, '{')
    comma = false
    for p in T.parameters
        s = string(p)
        if position(buf) + sizeof(s) > 20
            comma || print(buf, '…')
            break
        end
        comma && print(buf, ',')
        comma = true
        print(buf, s)
    end
    print(buf, '}')
    return String(take!(buf))
end

#Get typeMap for Jupyter completions
function complete_types(comps)
    typeMap = []
    for c in comps
        ctype = ""
        if !isempty(searchsorted(sorted_keywords, c))
            ctype = "keyword"
        elseif startswith(c, "\\:")
            ctype = get(emoji_symbols, c, "")
            isempty(ctype) || (ctype = "emoji: $ctype")
        elseif startswith(c, "\\")
            ctype = get(latex_symbols, c, "")
        else
            expr = Meta.parse(c, raise=false)
            if typeof(expr) == Symbol
                try
                    ctype = complete_type(Core.eval(current_module[], :(typeof($expr))))
                catch
                end
            elseif !isa(expr, Expr)
                ctype = complete_type(expr)
            elseif expr.head == :macrocall
                ctype = "macro"
            end
        end
        isempty(ctype) || push!(typeMap, Dict("text" => c, "type" => ctype))
    end
    return typeMap
end

function complete_request(socket, msg)
    code = msg.content["code"]
    cursor_chr = msg.content["cursor_pos"]
    cursorpos = chr2ind(msg, code, cursor_chr)
    if all(isspace, code[1:cursorpos])
        send_ipython(requests[], msg_reply(msg, "complete_reply",
                                 Dict("status" => "ok",
                                              "metadata" => Dict(),
                                              "matches" => String[],
                                              "cursor_start" => cursor_chr,
                                              "cursor_end" => cursor_chr)))
        return
    end

    codestart = find_parsestart(code, cursorpos)
    comps_, positions = REPLCompletions.completions(code[codestart:end], cursorpos-codestart+1, current_module[])
    comps = unique!(REPLCompletions.completion_text.(comps_)) # julia#26930
    # positions = positions .+ (codestart - 1) on Julia 0.7
    positions = (first(positions) + codestart - 1):(last(positions) + codestart - 1)
    metadata = Dict()
    if isempty(comps)
        # issue #530: REPLCompletions returns inconsistent results
        # for positions when no completions are found
        cursor_start = cursor_end = cursor_chr
    elseif isempty(positions) # true if comps to be inserted without replacement
        cursor_start = (cursor_end = ind2chr(msg, code, last(positions)))
    else
        cursor_start = ind2chr(msg, code, prevind(code, first(positions)))
        cursor_end = ind2chr(msg, code, last(positions))
        metadata["_jupyter_types_experimental"] = complete_types(comps)
    end
    send_ipython(requests[], msg_reply(msg, "complete_reply",
                                     Dict("status" => "ok",
                                                  "matches" => comps,
                                                  "metadata" => metadata,
                                                  "cursor_start" => cursor_start,
                                                  "cursor_end" => cursor_end)))
end

function kernel_info_request(socket, msg)
    send_ipython(requests[],
                 msg_reply(msg, "kernel_info_reply",
                           Dict("protocol_version" => "5.0",
                                        "implementation" => "ijulia",
                                        # TODO: "implementation_version" => IJulia version string from Pkg
                                        "language_info" =>
                                        Dict("name" => "julia",
                                             "version" =>
                                             string(VERSION.major, '.',
                                                    VERSION.minor, '.',
                                                    VERSION.patch),
                                             "mimetype" => "application/julia",
                                             "file_extension" => ".jl"),
                                        "banner" => "Julia: A fresh approach to technical computing.",
                                        "help_links" => [
                                                         Dict("text"=>"Julia Home Page",
                                                              "url"=>"http://julialang.org/"),
                                                         Dict("text"=>"Julia Documentation",
                                                              "url"=>"http://docs.julialang.org/"),
                                                         Dict("text"=>"Julia Packages",
                                                              "url"=>"http://pkg.julialang.org/")
                                                        ])))
end

function connect_request(socket, msg)
    send_ipython(requests[],
                 msg_reply(msg, "connect_reply",
                           Dict("shell_port" => profile["shell_port"],
                                        "iopub_port" => profile["iopub_port"],
                                        "stdin_port" => profile["stdin_port"],
                                        "hb_port" => profile["hb_port"])))
end

function shutdown_request(socket, msg)
    send_ipython(requests[], msg_reply(msg, "shutdown_reply",
                                     msg.content))
    sleep(0.1) # short delay (like in ipykernel), to hopefully ensure shutdown_reply is sent
    exit()
end

docdict(s::AbstractString) = display_dict(Core.eval(Main, helpmode(devnull, s)))

import Base: is_id_char, is_id_start_char
function get_token(code, pos)
    # given a string and a cursor position, find substring to request
    # help on by:
    #   1) searching backwards, skipping invalid identifier chars
    #        ... search forward for end of identifier
    #   2) search backwards to find the biggest identifier (including .)
    #   3) if nothing found, do return empty string
    # TODO: detect operators?

    startpos = pos
    while startpos > firstindex(code)
        if is_id_char(code[startpos])
            break
        else
            startpos = prevind(code, startpos)
        end
    end
    endpos = startpos
    while startpos >= firstindex(code) && (is_id_char(code[startpos]) || code[startpos] == '.')
        startpos = prevind(code, startpos)
    end
    startpos = startpos < pos ? nextind(code, startpos) : pos
    if !is_id_start_char(code[startpos])
        return ""
    end
    while endpos < lastindex(code) && is_id_char(code[endpos])
        endpos = nextind(code, endpos)
    end
    if !is_id_char(code[endpos])
        endpos = prevind(code, endpos)
    end
    return code[startpos:endpos]
end

function inspect_request(socket, msg)
    try
        code = msg.content["code"]
        s = get_token(code, chr2ind(msg, code, msg.content["cursor_pos"]))
        if isempty(s)
            content = Dict("status" => "ok", "found" => false)
        else
            d = docdict(s)
            content = Dict("status" => "ok",
                           "found" => !isempty(d),
                           "data" => d)
        end
        send_ipython(requests[], msg_reply(msg, "inspect_reply", content))
    catch e
        content = error_content(e, backtrace_top=:inspect_request);
        content["status"] = "error"
        send_ipython(requests[],
                     msg_reply(msg, "inspect_reply", content))
    end
end

function history_request(socket, msg)
    # we will just send back empty history for now, pending clarification
    # as requested in ipython/ipython#3806
    send_ipython(requests[],
                 msg_reply(msg, "history_reply",
                           Dict("history" => [])))
end

function is_complete_request(socket, msg)
    ex = Meta.parse(msg.content["code"], raise=false)
    status = Meta.isexpr(ex, :incomplete) ? "incomplete" : Meta.isexpr(ex, :error) ? "invalid" : "complete"
    send_ipython(requests[],
                 msg_reply(msg, "is_complete_reply",
                           Dict("status"=>status, "indent"=>"")))
end

const handlers = Dict{String,Function}(
    "execute_request" => execute_request,
    "complete_request" => complete_request,
    "kernel_info_request" => kernel_info_request,
    "inspect_request" => inspect_request,
    "connect_request" => connect_request,
    "shutdown_request" => shutdown_request,
    "history_request" => history_request,
    "is_complete_request" => is_complete_request,
    "comm_open" => comm_open,
    "comm_info_request" => comm_info_request,
    "comm_msg" => comm_msg,
    "comm_close" => comm_close
)


text:
using Distributed, MonteCarloObservable

const dqmc_file = joinpath(ENV["JULIA_DQMC"], "app/dqmc.jl")
const params_dir = joinpath(ENV["JULIA_DQMC"], "test/parameters")

function killworkers()
    while nprocs() > 1
        rmprocs(workers()[end])
    end
end

function change_sweeps_in_xml(xml, sweeps)
    str = read(xml, String)
    strs = split(str, "SWEEPS")
    start = findfirst(isequal('>'), strs[2])+1
    stop = findfirst(isequal('<'), strs[2])-1
    cur_sweeps = strs[2][start:stop]
    strs[2] = replace(strs[2], cur_sweeps => string(sweeps), count=1)
    write(xml, join(strs, "SWEEPS"))
    nothing
end

"""
Creates a fresh temporary directory for e.g. dqmc simulations.

To be used as

        intmpdir() do tmp_dir
            # whatever
        end
"""
function intmpdir(f)
    killworkers()
    mktempdir() do tmp_dir
        cd(tmp_dir) do

            # code
            tmp_dir = pwd()
            f(tmp_dir)

        end
    end
    killworkers()
end

"""
Run a dqmc simulation on a worker process. Returns the elapsed seconds.
"""
function run_dqmc(in_xml)
    addprocs(1) # Create worker to run dqmc.jl
    t = @elapsed @sync @spawnat workers()[1] begin
        empty!(ARGS)
        push!(ARGS, in_xml)
        include(dqmc_file)
    end
    killworkers()
    return t
end

"""
Temporarily set an environmental variable to a value
"""
function with_env_var(f, env, val)
    before_val = haskey(ENV, env) ? ENV[env] : nothing
    @info "Setting ENV[\"$env\"] = $val."
    ENV[env] = val
    f()
    if !isnothing(before_val)
        ENV[env] = before_val
    end
end

if !isfile(dqmc_file)
    @warn "dqmc.jl not found. Skipping Resume DQMC tests."
elseif !isdir(params_dir)
    @warn "parameters directory not found. Skipping Resume DQMC tests."
else
    @testset "Resume DQMC" begin

        intmpdir() do tmp_dir
            # cp(joinpath(params_dir, "resumelong.in.xml"), joinpath(tmp_dir, "resumelong.in.xml"))
            # cp(joinpath(params_dir, "resumeshort.in.xml"), joinpath(tmp_dir, "resumeshort.in.xml"))
            cp(joinpath(params_dir, "resumelong.in.xml"), "resumelong.in.xml")
            cp(joinpath(params_dir, "resumeshort.in.xml"), "resumeshort.in.xml")

            # Run resumelong
            println("Starting DQMC (resumelong)")
            tlong = run_dqmc("resumelong.in.xml")
            println("Done (resumelong)")

            # extract minutes and seconds from tlong/2
            # s = Second(floor(Int, tlong/2))
            # d = Dates.canonicalize(Dates.CompoundPeriod(s))
            # minutes = d.periods[1].value
            # seconds = d.periods[2].value

            @testset "Resume (extension run)" begin
                # Run resumeshort
                println("Starting DQMC (resumeshort)")
                run_dqmc("resumeshort.in.xml")
                println("Done (resumeshort)")

                mv("resumeshort.out.h5", "resumeshort.out.h5.running")
                change_sweeps_in_xml("resumeshort.in.xml", 200)

                # Run resumeshort once again
                println("Starting DQMC (resumeshort - again)")
                run_dqmc("resumeshort.in.xml")
                println("Done (resumeshort - again)")

                # Comparison
                clong = ts_flat("resumelong.out.h5", "obs/configurations");
                cshort = ts_flat("resumeshort.out.h5", "obs/configurations");
                l = Int(size(clong, ndims(clong)) / 2)
                @test @views isapprox(clong[:,:,:,1:l], cshort[:,:,:,1:l]) # should trivially pass
                @test isapprox(clong, cshort)

                glong = ts_flat("resumelong.out.h5", "obs/greens");
                gshort = ts_flat("resumeshort.out.h5", "obs/greens");
                @test @views isapprox(glong[:,:,1:l], gshort[:,:,1:l]) # should trivially pass
                @test isapprox(glong, gshort)
            end

            # VERY ERROR PRONE....
            # @testset "Resume after safe exit (WTL)" begin
            #     cp("resumelong.in.xml", "resumewtl.in.xml")

            #     with_env_var("WALLTIMELIMIT", "0-00:00:$(floor(Int, tlong/2))") do
            #         try
            #             println("Starting DQMC (resumetwtl 1)")
            #             run_dqmc("resumewtl.in.xml")
            #             println("Done (resumetwtl 1)")
            #         catch err
            #             if !isa(err, ProcessExitedException)
            #                 throw(err)
            #             end
            #         end

            #         try
            #             println("Starting DQMC (resumetwtl 2)")
            #             run_dqmc("resumewtl.in.xml")
            #             println("Done (resumetwtl 2)")
            #         catch err
            #             if !isa(err, ProcessExitedException)
            #                 throw(err)
            #             end
            #         end

            #         # check that job hasn't finished but did dump something
            #         if isfile("resumewtl.out.h5")
            #             @warn "RESUME WTL: WTL too long, job did finish. Skipping."
            #         elseif !isfile("resumewtl.out.h5.running")
            #             @warn "RESUME WTL: Job didn't do anything? Skipping."
            #         elseif !("configurations" in listobs("resumewtl.out.h5.running"))
            #             @warn "RESUME WTL: Found .running file but no confs. WTL too short. Skipping."
            #         else
            #             c = loadobs_frommemory("resumewtl.out.h5.running", "obs/configurations")
            #             clong = loadobs_frommemory("resumelong.out.h5", "obs/configurations")
            #             @test isapprox(c, clong[:,:,:,1:length(c)])
            #         end
            #     end
            # end
        end # intmpdir

    end # resume testset
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
mutable struct DLUE{V<:Real}
    F::Matrix{V}
    B::Matrix{V}
    C::OutputMatrix{V}
    L::Matrix{V}
    E::Matrix{V}
    z::Vector{V}

    function DLUE(A::T, B::T, C::OutputMatrix{V}, L::T) where {V<:Real, T<:Matrix{V}}
        new{V}(A - L*C,
            B,
            C,
            L,
            zeros(size(A,1), 0),
            zeros(V, size(A,1)))
    end
end

function DLUE(s::LdtiSubsystem{T}, L::AbstractMatrix{T}) where {T}
    o = DLUE(s.A, s.B, s.C, Matrix(L))
    setNeighbourCoupling(o, s.E)
    o
end

function setInitialState(o::DLUE{V}, z0::Vector{V}) where {V<:Real}
    setfield!(o, :z, z0)
end

function setNeighbourCoupling(o::DLUE{V}, coupling::AbstractMatrix{V}) where {V}
    setfield!(o, :E, isempty(coupling) ? Matrix{V}(undef, o.nx, 0) : coupling)
end

function updateState(o::DLUE{V}, u::Vector{V}, y::Vector{V}, hatbarx::AbstractVector{V}) where {V<:Real}
    if isempty(hatbarx)
        setfield!(o, :z, o.F*o.z + o.B*u + o.L*y)
    else
        setfield!(o, :z, o.F*o.z + o.B*u + o.L*y + o.E*hatbarx)
    end
end

updateState(o::DLUE{V}, u::V, y::Vector{V}, hatbarx::Vector{V}) where {V} = updateState(o, [u], y, hatbarx)

outputMap(o::DLUE) = o.C*o.z

function Base.getproperty(o::DLUE, s::Symbol)
    if s === :nx
        return size(o.F, 1)
    else
        getfield(o, s)
    end
end

function Base.setproperty!(::DLUE, ::Symbol, val)
    error("private object")
end


text:
using RecipesBase, Unitful
import Measurements: value

@testset "Plotting" begin
    wave = range(1e4, 3e4, length = 1000)
    sigma = 0.1 .* sin.(wave)
    T = 6700
    flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ± sigma
    wave *= u"angstrom"
    flux *= u"erg/s/cm^2/angstrom"
    spec = spectrum(wave, flux, name = "Test Spectrum")

    rec = RecipesBase.apply_recipe(Dict{Symbol,Any}(), spec)
    @test getfield(rec[1], 1) == Dict{Symbol,Any}(:label => "",
        :xlabel => "wave",
        :ylabel => "flux density",
        :seriestype => :step)
    @test rec[1].args == (ustrip.(spec.wave), value.(ustrip.(spec.flux)))

    strip_spec = ustrip(spec)

    rec = RecipesBase.apply_recipe(Dict{Symbol,Any}(), strip_spec)
    @test getfield(rec[1], 1) == Dict{Symbol,Any}(:label => "",
        :xlabel => "wave",
        :ylabel => "flux density",
        :seriestype => :step)
    @test rec[1].args == (strip_spec.wave, value.(strip_spec.flux))
end


text:
module Covering

import Voronoi

abstract type Element end

mutable struct Edge <: Element
   origin::Tuple{Real, Real}
   next::Union{Element, Nothing}
end

mutable struct Arc <: Element
   origin::Tuple{Real, Real}
   next::Union{Element, Nothing}
end

mutable struct Section
   center::Union{Tuple{Real, Real}, Nothing}
   borderHead::Union{Element, Nothing}
end

function Section(points::Array{Tuple{Real, Real}, 1})
    n = size(points)[1]

    edge = head = Edge(points[1], nothing)

    for i in 2:n
        edge.next = Edge(points[i], nothing)
        edge = edge.next
    end

    edge.next = head

    return Section((0, 0), head)
end

function Section(section::Section)
    edge = head = Edge(section.borderHead.origin, nothing)

    copying = section.borderHead.next
    while copying != section.borderHead
        edge.next = Edge(copying.origin, nothing)
        edge = edge.next
        copying = copying.next
    end

    edge.next = head

    return Section((0, 0), head)
end

function segmentArcIntersections(segment::Voronoi.Diagram.HalfEdge, center::Tuple{Real, Real}, r::Real)
    return segmentArcIntersections(segment.origin, segment.next.origin, center, r)
end

function segmentArcIntersections(p::Tuple{Real, Real}, q::Tuple{Real, Real}, center::Tuple{Real, Real}, r::Real)
    function dot(a, b)
        return a[1]*b[1] + a[2]*b[2]
    end

    PQ = Voronoi.Geometry.subVector(q, p)
    CP = Voronoi.Geometry.subVector(p, center)

    a = dot(PQ, PQ)
    b = 2*dot(CP, PQ)
    c = dot(CP, CP) - r^2

    Δ = b^2 - 4a*c # discriminant

    if Δ < 0 # no intersection
        return []
    elseif Δ == 0 # one intersection
        α = (-b - sqrt(Δ))/2a

        if 0 <= α <= 1
            return [(p[1] + α*PQ[1], p[2] + α*PQ[2])]
        else
            return []
        end
    else # two intersections
        α1 = (-b - sqrt(Δ))/2a
        α2 = (-b + sqrt(Δ))/2a

        if (α1 < α2)
            if 0 <= α1 <= α2 <= 1
                return [(p[1] + α1*PQ[1], p[2] + α1*PQ[2]), (p[1] + α2*PQ[1], p[2] + α2*PQ[2])]
            elseif 0 <= α1 <= 1
                return [(p[1] + α1*PQ[1], p[2] + α1*PQ[2])]
            elseif 0 <= α2 <= 1
                return [(p[1] + α2*PQ[1], p[2] + α2*PQ[2])]
            else
                return []
            end
        else
            if 0 <= α2 <= α1 <= 1
                return [(p[1] + α2*PQ[1], p[2] + α2*PQ[2]), (p[1] + α1*PQ[1], p[2] + α1*PQ[2])]
            elseif 0 <= α1 <= 1
                return [(p[1] + α1*PQ[1], p[2] + α1*PQ[2])]
            elseif 0 <= α2 <= 1
                return [(p[1] + α2*PQ[1], p[2] + α2*PQ[2])]
            else
                return []
            end
        end
    end
end

function distanceOutsideCircle(a, center, r)
    d = Voronoi.Geometry.distanceSquared(a, center)
    return d - r^2
end

function intersectedSegment(segment::Voronoi.Diagram.HalfEdge, center::Tuple{Real, Real}, r::Real)
    p = segment.origin
    q = segment.next.origin

    docp = distanceOutsideCircle(p, center, r)
    docq = distanceOutsideCircle(q, center, r)

    if docp <= 0 && docq <= 0
        el = head = Edge(p, nothing)
    else
        intersections = segmentArcIntersections(segment, center, r)

        if docq >= 0 && docp >= 0 && length(intersections) == 1
            el = head = Arc(p, nothing)
        else
            if docp < 0
                inside = true
                el = head = Edge(p, nothing)
            else
                inside = false
                el = head = Arc(p, nothing)
            end

            for i in intersections
                if inside
                    el.next = Arc(i, nothing);
                else
                    el.next = Edge(i, nothing);
                end

                el = el.next
                inside = !inside
            end
        end
    end

    return head, el
end

function voronoiRegionToCoveringSection(region::Voronoi.Diagram.Region, r::Real)
    he = region.borderHead.next
    head, tail = intersectedSegment(he, region.generator, r)

    el = head
    while (he != region.borderHead)
        he = he.next
        tail.next, tail = intersectedSegment(he, region.generator, r)
    end
    tail.next = head

    return Section(region.generator, head)
end

mutable struct Partition
    sections::Array{Covering.Section}
    r::Real
end

function voronoiDiagramToPartition(V::Voronoi.Diagram.DCEL, r::Real)
    n = length(V.regions)
    sections = Array{Section, 1}(undef, n)

    for i in 1:n
        sections[i] = voronoiRegionToCoveringSection(V.regions[i], r)
    end

    return Partition(sections, r)
end

function radianToDegrees(rad)
    return rad*180/π
end

function segmentAngle(p, q)
    num = q[2] - p[2]
    den = q[1] - p[1]
    if num == 0 && den > 0
        angle = 0
    elseif num == 0 && den < 0
        angle = π
    elseif den == 0 && num > 0
        angle = π/2
    elseif den == 0 && num < 0
        angle = 3π/2
    elseif num > 0 && den > 0
        angle = atan(num/den)
    elseif num > 0 && den < 0
        angle = π - atan(-num/den)
    elseif num < 0 && den < 0
        angle = π + atan(num/den)
    else # num < 0 && den > 0
        angle = -atan(-num/den)
    end

    return angle
end

function areaAndGradient(section::Section, r::Real)
    if section.borderHead == nothing
        return 0, 0, (0, 0)
    end

    coveredArea = 0
    el = section.borderHead
    gᵣ = 0
    gₛx = gₛy = 0.0

    firstTime = true # gimme do whiles, Julia, dammit!
    while (el != section.borderHead || firstTime)
        firstTime = false
        if (el isa Arc)
            θ1 = segmentAngle(section.center, el.origin)
            θ2 = segmentAngle(section.center, el.next.origin)
            θ = θ2 - θ1
            if θ < 0
                θ += 2π
            end

            coveredArea += section.center[1]*r*(sin(θ2) - sin(θ1))
            coveredArea += ((r^2)/2)*(θ + sin(θ2)*cos(θ2) - sin(θ1)*cos(θ1))

            gᵣ += -r*θ # (θ/(2π))*2πr
            gₛx += r*(sin(θ1) - sin(θ2))
            gₛy += r*(cos(θ2) - cos(θ1))
        else # isa Edge
            v = el.origin
            w = el.next.origin

            coveredArea += (w[2] - v[2])*(v[1] + w[1])/2
        end

        el = el.next
    end

    return coveredArea, gᵣ, (gₛx, gₛy)
end

function areaAndGradient(P::Partition)
    n = length(P.sections)

    coveredArea = gᵣ = 0
    gₛ = fill((0.0, 0.0), n)

    for i in 1:n
        a, g_r, g_s_i = areaAndGradient(P.sections[i], P.r)
        coveredArea += a
        gᵣ += g_r
        gₛ[i] = g_s_i
    end

    return coveredArea, gᵣ, gₛ
end

end # module

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
"""
    embedding(small::AbstractMesh, big::AbstractMesh) -> S::AbstractSparseArray
    predicate(a,b) => extra constraints that must both be fullfiled by each
    cell `a` and `b` that are found to exist in `small` and `big` respectively.
    `a` and `b` are the indices of the corresponding edges. It is the user's
    responsibility to ensure that the predicate allows for all cells of `small`
    to be mapped to a cell in `big`. Specifying a predicate is useful for cases
    where cells in `small` matches/overlaps with multiple cells in `big`
    #TODO: change pred to accept simplices
"""
function embedding(small::AbstractMesh, big::AbstractMesh; predicate=(a,b)->true)

    @assert dimension(small) == 1
    @assert dimension(big)   == 1

    @assert length(big) > 0
    @assert length(small) > 0

    T = coordtype(small)
    tol = sqrt(eps(T))

    # store the edges of the big mesh in an octree
    ctrs_big = [cartesian(center(chart(big,v))) for v in big]
    tree_big = Octree(ctrs_big, zeros(size(ctrs_big)))

    # rows = 1:length(small)
    cols = zeros(length(small))
    sgns = zeros(length(small))

    cells_big = collect(cells(big))
    for (i,c) in enumerate(small)
        V = chart(small, c).vertices
        ch = chart(small,c)
        x = cartesian(center(ch))
        pred(c,s) = fitsinbox(x, 0.0, c, s+tol)
        found = false
        for box in boxes(tree_big, pred)
            for j in box
                c′ = cells_big[j]
                V′ = chart(big,c′).vertices
                ch′ = chart(big, c′)
                if (norm(V[1]-V′[1]) + norm(V[2]-V′[2]) < tol) && predicate(i,j)
                    cols[i] = j
                    sgns[i] = +1
                    found = true
                    break
                elseif (norm(V[1]-V′[2]) + norm(V[2]-V′[1]) < tol) && predicate(i,j)
                    cols[i] = j
                    sgns[i] = -1
                    found = true
                    break
                end
            end
            found && break
        end
    end

    @assert !any(cols .== 0) "Unmapped cells exist"
    @assert !any(sgns .== 0)
    sparse(1:length(small), cols, sgns, length(small), length(big))
end

function embedding_topo(small::AbstractMesh, big::AbstractMesh)

    @assert dimension(small) == dimension(big) == 1
    @assert vertices(big) == vertices(small)

    @assert length(big) * length(small) > 0

    T = coordtype(small)
    tol = sqrt(eps(T))

    # store the edges of the big mesh in an octree
    ctrs_big = [cartesian(center(chart(big,v))) for v in big]
    tree_big = Octree(ctrs_big, zeros(size(ctrs_big)))

    # rows = 1:length(small)
    cols = zeros(length(small))
    sgns = zeros(length(small))

    cells_big = collect(cells(big))
    for (i,c) in enumerate(small)
        V = chart(small, c).vertices
        ch = chart(small,c)
        x = cartesian(center(ch))
        pred(c,s) = fitsinbox(x, 0.0, c, s+tol)
        found = false
        for box in boxes(tree_big, pred)
            for j in box
                c′ = cells_big[j]

                if c′ == c
                    cols[i] = j
                    sgns[i] = +1
                    found = true
                    break
                elseif c′ == reverse(c)
                    cols[i] = j
                    sgns[i] = -1
                    found = true
                    break
                end

                # V′ = chart(big,c′).vertices
                # ch′ = chart(big, c′)
                # if (norm(V[1]-V′[1]) + norm(V[2]-V′[2]) < tol) && predicate(i,j)
                #     cols[i] = j
                #     sgns[i] = +1
                #     found = true
                #     break
                # elseif (norm(V[1]-V′[2]) + norm(V[2]-V′[1]) < tol) && predicate(i,j)
                #     cols[i] = j
                #     sgns[i] = -1
                #     found = true
                #     break
                # end

            end
            found && break
        end
    end

    @assert !any(cols .== 0) "Unmapped cells exist"
    @assert !any(sgns .== 0)
    sparse(1:length(small), cols, sgns, length(small), length(big))
end


text:
using GPUifyLoops

function kernel(A)
    @loop for i in (1:size(A,1);
                    threadIdx().x)
        A[i] = 2*A[i]
    end
    @synchronize
end

data = Array{Float32}(undef, 1024)
kernel(data)

@static if Base.find_package("CuArrays") !== nothing
    using CuArrays
    using CUDAnative

    kernel(A::CuArray) = @launch CUDA() kernel(A, threads=length(A))

    data = CuArray{Float32}(undef, 1024)
    kernel(data)
end


text:
function _has_lb(model::Optimizer, index::MOI.VariableIndex)
    return _info(model, index).lb !== nothing
end
function _has_ub(model::Optimizer, index::MOI.VariableIndex)
    return _info(model, index).ub !== nothing
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) && _has_ub(model, index)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) && _has_lb(model, index)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.Interval{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) &&
           _has_lb(model, index) &&
           _has_ub(model, index)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) &&
           _info(model, index).lb == _info(model, index).ub
end

function _build_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.EqualTo{T},
) where {T <: Integer}
    v = _info(model, f).variable
    return XeqC(v, s.value)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne},
)
    index = MOI.VariableIndex(c.value)
    if !MOI.is_valid(model, index)
        return false
    end

    return _info(model, index).type == BINARY
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.Integer},
)
    index = MOI.VariableIndex(c.value)
    if !MOI.is_valid(model, index)
        return false
    end

    return _info(model, index).type == INTEGER
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{MOI.VariableIndex, <:Any},
)
    MOI.throw_if_not_valid(model, c)
    return MOI.VariableIndex(MOI.VariableIndex(c.value))
end

function MOI.set(
    ::Optimizer,
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.VariableIndex, S},
    ::MOI.VariableIndex,
) where {S}
    throw(MOI.SettingVariableIndexFunctionNotAllowed())
    return
end

# CP.Domain.
function MOI.supports_constraint(
    ::Optimizer,
    ::MOI.VariableIndex,
    ::CP.Domain{T},
) where {T <: Integer}
    return true
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::CP.Domain{T},
) where {T <: Integer}
    # Add each possible value separately. Same logic as
    # https://github.com/radsz/jacop/blob/develop/src/main/java/org/jacop/fz/VariablesParameters.java#L238-L267
    v = _info(model, f).variable
    for e in s.values
        jcall(v, "addDom", (jint, jint), e, e)
    end

    index = MOI.ConstraintIndex{F, S}(length(model.constraint_info) + 1)
    model.constraint_info[index] = ConstraintInfo(index, nothing, f, s)
    return index
end


text:
stress_scale = 1.0e+5
strain_scale = 1

# tid = parse(Int64, ARGS[1])
force_scale = 2.0
tid=200

printstyled("tid=$tid\n", color=:green)

testtype = "NeuralNetwork2D"
nntype = "stiffmat"
include("nnutil.jl")
printstyled("idx = $idx, tid=$tid\n", color=:green)

H0 = [1.04167e6  2.08333e5  0.0
      2.08333e5  1.04167e6  0.0
      0.0        0.0        4.16667e5]/stress_scale

#s = ae_to_code("Data/NNLear.mat", nntype)
#s = ae_to_code("Data/NNPreLSfit_$(idx).mat", nntype)
i = 5
#nnname="Data/$(nntype)/NN_Train_$(idx)_ite$(i).mat"
#nnname="Data/$(nntype)/NNLearn_$(idx)_ite$(i).mat"
nnname="Data/$(nntype)/NNPreLSfit_$(idx)_spd_Chol_Orth_$(i).mat"
s = ae_to_code(nnname, nntype)

eval(Meta.parse(s))
# density 4.5*(1 - 0.25) + 3.2*0.25
#fiber_fraction = 0.25
#todo
prop = Dict("name"=> testtype, "rho"=> 0.8000, "nn"=>post_nn)

T = 0.2
NT = 200

# nx_f, ny_f = 12, 4
# homogenized computaional domain
# number of elements in each directions
nx, ny = 10, 5
Lx, Ly = 0.1, 0.05 #m
porder = 2

nodes, EBC, g, gt, FBC, fext, ft = BoundaryCondition(tid, nx, ny,porder, Lx, Ly;force_scale=force_scale)

ndofs=2
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, 3))
    end
end

domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

ts = LinRange(0, T, NT+1)
adaptive_solver_args = Dict("Newmark_rho"=> 0.0,
                          "Newton_maxiter"=>10,
                          "Newton_Abs_Err"=>1e-4,
                          "Newton_Rel_Err"=>1e-6,
                          "damped_Newton_eta" => 1.0)

globdat, domain, ts = AdaptiveSolver("NewmarkSolver", globdat, domain, T, NT, adaptive_solver_args)

close("all")
visσ(domain)
axis("equal")
savefig("Debug/order$porder/test_stress$(tid)_$(idx)_$force_scale.png")

close("all")
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], ny*porder+1, nx*porder+1)[end,end] for i = 1:length(domain.history["state"])]
plot(ts, ux)
savefig("Debug/order$porder/test_ux$(tid)_$(idx)_$force_scale.png")

close("all")
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], ny*porder+1, nx*porder+1)[end,end] for i = 1:length(domain.history["state"])]
plot(ts, uy)
savefig("Debug/order$porder/test_uy$(tid)_$(idx)_$force_scale.png")