text:
#=
    - Julia Version : 1.4.0
    - Author        : Alexanders Bills and Abhinav Misalkar
    - Date          : 14th May 2020
    - Title         : Overhead Costs and Summary
=#

function OEM(input_OEM,cost,breakdown,materials_breakdown)

    cost_per_cell_positive_active_material              = input_OEM[1]                       ##   $/cell
    cost_per_cell_negative_active_material              = input_OEM[2]                       ##   $/cell
    cost_per_cell_positive_conductive                   = input_OEM[3]                       ##   $/cell
    cost_per_cell_positive_binder                       = input_OEM[4]                       ##   $/cell
    cost_per_cell_binder_solvent_NMP                    = input_OEM[5]                       ##   $/cell
    cost_per_cell_negative_conductive                   = input_OEM[6]                       ##   $/cell
    cost_per_cell_negative_binder                       = input_OEM[7]                       ##   $/cell
    cost_per_cell_positive_current_collector            = input_OEM[8]                       ##   $/cell
    cost_per_cell_negative_current_collector            = input_OEM[9]                       ##   $/cell
    cost_per_cell_separator                             = input_OEM[10]                      ##   $/cell
    cost_per_cell_electrolyte                           = input_OEM[11]                      ##   $/cell
    cost_per_cell_positive_terminal_assembly            = input_OEM[12]                      ##   $/unit
    cost_per_cell_negative_terminal_assembly            = input_OEM[13]                      ##   $/unit
    cost_per_cell_cell_container                        = input_OEM[14]                      ##   $/unit

    receiving_direct_labor                              = input_OEM[15]                      ##   hours/year
    positive_material_preparation_direct_labor          = input_OEM[16]                      ##   hours/year
    negative_material_preparation_direct_labor          = input_OEM[17]                      ##   hours/year
    positive_electrode_coating_direct_labor             = input_OEM[18]                      ##   hours/year
    negative_electrode_coating_direct_labor             = input_OEM[19]                      ##   hours/year
    binder_solvent_NMP_recovery_direct_labor            = input_OEM[20]                      ##   hours/year
    positive_calendering_direct_labor                   = input_OEM[21]                      ##   hours/year
    negative_calendering_direct_labor                   = input_OEM[22]                      ##   hours/year
    inter_process_material_handling_direct_labor        = input_OEM[23]                      ##   hours/year
    electrode_slitting_direct_labor                     = input_OEM[24]                      ##   hours/year
    vacuum_drying_direct_labor                          = input_OEM[25]                      ##   hours/year
    control_laboratory_direct_labor                     = input_OEM[26]                      ##   hours/year
    cell_stacking_direct_labor                          = input_OEM[27]                      ##   hours/year
    current_collector_welding_direct_labor              = input_OEM[28]                      ##   hours/year
    cell_insertion_container_direct_labor               = input_OEM[29]                      ##   hours/year
    electrolyte_filling_and_cell_sealing_direct_labor   = input_OEM[30]                      ##   hours/year
    dry_room_control_direct_labor                       = input_OEM[31]                      ##   hours/year
    formation_cycling_direct_labor                      = input_OEM[32]                      ##   hours/year
    cell_sealing_rack_unloading_direct_labor            = input_OEM[33]                      ##   hours/year
    charge_retention_testing_direct_labor               = input_OEM[34]                      ##   hours/year
    rejected_cell_scrap_recycle_assembly_direct_labor   = input_OEM[35]                      ##   hours/year
    shipping_direct_labor                               = input_OEM[36]                      ##   hours/year

    total_plant_area                                    = input_OEM[37]                      ##   m2
    total_capital_equipment                             = input_OEM[38]                      ##   USD Million
    total_units_mfg_per_year                            = input_OEM[39]                      ##   units/year
    cell_energy_storage                                 = input_OEM[40]                      ##   kWh
    cost_verbosity                                      = input_OEM[41]                      ##   1 or 0
    mfg_capacity                                        = input_OEM[42]                      ## MWh/year
    system                                              = input_OEM[43]                      ## Cell or Pack
    if system == 1
        system = "Cell"
    elseif system == 2
        system = "Pack"
    end

    ## Investment Costs
    building_investment = (cost.general_costs.building * total_plant_area)

                            ###################### Cost Breakdown Analysis #####################

    total_cost_carbon_and_binders =
            cost_per_cell_positive_conductive  +
            cost_per_cell_positive_binder      +
            cost_per_cell_binder_solvent_NMP   +
            cost_per_cell_negative_conductive  +
            cost_per_cell_negative_binder

    total_cost_materials =
            cost_per_cell_positive_active_material      +
            cost_per_cell_negative_active_material      +
            total_cost_carbon_and_binders               +
            cost_per_cell_positive_current_collector    +
            cost_per_cell_negative_current_collector    +
            cost_per_cell_separator                     +
            cost_per_cell_electrolyte

    total_purchased_items =
            cost_per_cell_positive_terminal_assembly    +
            cost_per_cell_negative_terminal_assembly    +
            cost_per_cell_cell_container

    total_unit_cost_all_materials =
            total_cost_materials + total_purchased_items

                                            ### Direct Labor Cost Breakdown ###

    direct_labor_per_unit_electrode_processing =
        (cost.general_costs.dir_labor / total_units_mfg_per_year) * (
            positive_material_preparation_direct_labor +
            negative_material_preparation_direct_labor +
            positive_electrode_coating_direct_labor +
            negative_electrode_coating_direct_labor +
            binder_solvent_NMP_recovery_direct_labor +
            positive_calendering_direct_labor +
            negative_calendering_direct_labor +
            inter_process_material_handling_direct_labor +
            electrode_slitting_direct_labor +
            vacuum_drying_direct_labor)

    direct_labor_per_unit_cell_assembly =
        (cost.general_costs.dir_labor / total_units_mfg_per_year) * (
            cell_stacking_direct_labor +
            current_collector_welding_direct_labor +
            cell_insertion_container_direct_labor +
            electrolyte_filling_and_cell_sealing_direct_labor +
            dry_room_control_direct_labor)

    direct_labor_per_unit_formation_cycling_testing_sealing =
        (cost.general_costs.dir_labor / total_units_mfg_per_year) * (
            formation_cycling_direct_labor +
            charge_retention_testing_direct_labor +
            cell_sealing_rack_unloading_direct_labor)

    direct_labor_per_unit_cell_materials_rejection_recycling =
        (cost.general_costs.dir_labor / total_units_mfg_per_year) *
        (rejected_cell_scrap_recycle_assembly_direct_labor)

    direct_labor_per_unit_receiving_shipping =
        (cost.general_costs.dir_labor / total_units_mfg_per_year) *
        (receiving_direct_labor + shipping_direct_labor)

    direct_labor_per_unit_control_laboratory =
        (cost.general_costs.dir_labor / total_units_mfg_per_year) *
        (control_laboratory_direct_labor)

    total_direct_labor_per_unit =
        direct_labor_per_unit_electrode_processing                  +
        direct_labor_per_unit_cell_assembly                         +
        direct_labor_per_unit_formation_cycling_testing_sealing     +
        direct_labor_per_unit_cell_materials_rejection_recycling    +
        direct_labor_per_unit_receiving_shipping                    +
        direct_labor_per_unit_control_laboratory

                                    ###################### Overhead Costs #####################
    depreciation_per_unit =
        ((cost.ovhd_rate.depr_capital_equip * total_capital_equipment +                 # total_capital_equipment currently USD
          cost.ovhd_rate.depr_building_inv * building_investment)                       # building_investment currently USD
        ) / total_units_mfg_per_year

    variable_overhead_per_unit =
        total_direct_labor_per_unit * cost.ovhd_rate.var_ovhd_dir_labor +
        depreciation_per_unit * cost.ovhd_rate.var_ovhd_depr

    GSA_per_unit =
        cost.ovhd_rate.GSA_dir_labor_var_ovhd *
        (total_direct_labor_per_unit + variable_overhead_per_unit) +
        cost.ovhd_rate.GSA_depr * depreciation_per_unit

    R_and_D_per_unit = depreciation_per_unit * cost.ovhd_rate.RnD_depr

    total_fixed_expenses_per_unit =
        GSA_per_unit + R_and_D_per_unit + depreciation_per_unit

    total_variable_cost_per_unit =
        total_direct_labor_per_unit +
        variable_overhead_per_unit +
        total_unit_cost_all_materials

                                    ###################### Investment Costs #####################
    total_launch_cost =
        (
            cost.ovhd_rate.launch_dir_annual_materials                  *
            total_unit_cost_all_materials                               +
            cost.ovhd_rate.launch_dir_labor_var_ovhd                    *
            (total_direct_labor_per_unit + variable_overhead_per_unit)) *
            (total_units_mfg_per_year)

    # print("\ncost.ovhd_rate.launch_dir_annual_materials = ", cost.ovhd_rate.launch_dir_annual_materials)
    # print("\ntotal_unit_cost_all_materials = ", total_unit_cost_all_materials)
    # print("\ncost.ovhd_rate.launch_dir_labor_var_ovhd = ", cost.ovhd_rate.launch_dir_labor_var_ovhd)
    # print("\ntotal_direct_labor_per_unit = ", total_direct_labor_per_unit)
    # print("\nvariable_overhead_per_unit = ", variable_overhead_per_unit)
    # print("\ntotal_units_mfg_per_year = ", total_units_mfg_per_year)
    # print("\ntotal_launch_cost = ", total_launch_cost)

    total_working_capital =
        (total_variable_cost_per_unit    *
            total_units_mfg_per_year     *
            cost.ovhd_rate.working_capital_rate)

    total_investment =                                                                                # currently USD
        building_investment     +
        total_capital_equipment +                                                                       # currently USD
        total_launch_cost       +                                                                       # currently USD
        total_working_capital                                                                           # currently USD

    profits_after_taxes_per_unit =
        (cost.ovhd_rate.profit_after_taxes * total_investment) /                                        # USD
        total_units_mfg_per_year

    warrant_cost_per_unit =
        cost.ovhd_rate.warranty            * (
            total_unit_cost_all_materials  +
            total_direct_labor_per_unit    +
            variable_overhead_per_unit     +
            GSA_per_unit                   +
            R_and_D_per_unit               +
            depreciation_per_unit          +
            profits_after_taxes_per_unit)

                                        ###################### Cost to OEM #####################
    cell_price_to_OEM =
    total_unit_cost_all_materials                 +
    total_direct_labor_per_unit                   +
    variable_overhead_per_unit                    +
    GSA_per_unit                                  +
    R_and_D_per_unit                              +
    depreciation_per_unit                         +
    profits_after_taxes_per_unit                  +
    warrant_cost_per_unit

    dollars_per_kWh = cell_price_to_OEM / cell_energy_storage

    # # if cost_verbosity == 0
    # print("\n")
    # # println("\n\n************************************************ Results ************************************************")
    # print("\nNo of cells mfg (Million)  =  ", total_units_mfg_per_year)
    # print("\nProduction Capacity (MWh)  =  ", mfg_capacity)
    # print("\nJulia dollars/kWh          =  ", dollars_per_kWh)
    # print("\ntotal_cell_cost_to_OEM     =  ", cell_price_to_OEM)
    # end

    if cost_verbosity == 1
        print("\n\n************************************************ Invsetment Costs (Million USD) ************************************************ ")

        print("\n\nbuilding_investment       = ", building_investment)
        print("\ntotal_launch_cost           = ", total_launch_cost)
        print("\ntotal_working_capital       = ", total_working_capital)
        print("\ntotal_investment            = ", total_investment)

        print("\n\n************************************************ Unit Cost Battery Pack (USD) ************************************************ ")
        if system == "Pack"
            print("\nVariable Cost")
            print("\n\tMaterials and Purchased Items")
            print("\n\t\ttotal_cost_pack_cell_materials     = ",total_cost_pack_cell_materials)
            print("\n\t\ttotal_pack_purchased_items         = ",total_pack_purchased_items)
            print("\n\t\ttotal_cost_module_material         = ",total_cost_module_materials)
            print("\n\t\ttotal_cost_pack_materials          = ",total_cost_pack_materials)
            print("\n\t\ttotal_pack_cost_all_materials      = ",total_pack_cost_all_materials)
        else
            print("\nVariable Cost")
            print("\n\tMaterials and Purchased Items")
            print("\n\t\ttotal_cost_materials               = ",total_cost_materials)
            print("\n\t\ttotal_purchased_items              = ",total_purchased_items)
            print("\n\t\ttotal_unit_cost_all_materials      = ",total_unit_cost_all_materials)

        end

        print("\n\nDirect Labor @ USD", cost.general_costs.dir_labor , "/ hour")
        print("\n\t\tdirect_labor_per_unit_electrode_processing                = ",direct_labor_per_unit_electrode_processing , " ")
        print("\n\t\tdirect_labor_per_unit_cell_assembly                       = ",direct_labor_per_unit_cell_assembly , " ")
        print("\n\t\tdirect_labor_per_unit_formation_cycling_testing_sealing   = ",direct_labor_per_unit_formation_cycling_testing_sealing , " ")

        if system == "Pack"
            print("\n\t\tdirect_labor_per_unit_module_and_battery_assembly     = ",direct_labor_per_unit_module_and_battery_assembly , " ")
        end

        print("\n\t\tdirect_labor_per_unit_cell_materials_rejection_recycling  = ",direct_labor_per_unit_cell_materials_rejection_recycling , " ")
        print("\n\t\tdirect_labor_per_unit_receiving_shipping                  = ",direct_labor_per_unit_receiving_shipping , " ")
        print("\n\t\tdirect_labor_per_unit_control_laboratory                  = ",direct_labor_per_unit_control_laboratory , " ")
        print("\n\t\ttotal_direct_labor_per_unit                               = ",total_direct_labor_per_unit , " ")

        print("\n\nvariable_overhead_per_unit                                  = ",variable_overhead_per_unit , " ")
        print("\ntotal_variable_cost_per_unit                                  = ",total_variable_cost_per_unit , " ")

        print("\n\nFixed Expenses")
        print("\n\tGSA_per_unit                     = ",GSA_per_unit ," ")
        print("\n\tR_and_D_per_unit                 = ",R_and_D_per_unit ," ")
        print("\n\tdepreciation_per_unit            = ",depreciation_per_unit ," ")
        print("\n\ttotal_fixed_expenses_per_unit    = ",total_fixed_expenses_per_unit ," ")
        print("\n\tprofits_after_taxes_per_unit     = ",profits_after_taxes_per_unit ," ")

        print("\n\n************************************************ Summary of Unit Costs (USD) ************************************************ ")
        print("\ntotal_cost_pack_cell_materials     = ",total_cost_materials)

        if system == "Pack"
            print("\ntotal purchased items              = ", (cost_per_unit_hardware + total_cost_module_materials + total_cost_pack_materials) ," ")
        else
            print("\ntotal purchased items              = ", (total_purchased_items))
        end

        print("\ntotal_direct_labor_per_unit        = ",total_direct_labor_per_unit , " ")
        print("\nvariable_overhead_per_unit         = ",variable_overhead_per_unit , " ")
        print("\nGSA_per_unit                       = ",GSA_per_unit ," ")
        print("\nR_and_D_per_unit                   = ",R_and_D_per_unit ," ")
        print("\ndepreciation_per_unit              = ",depreciation_per_unit ," ")
        print("\nprofits_after_taxes_per_unit       =  ",profits_after_taxes_per_unit ," ")
        print("\nwarrant_cost_per_unit              = ",warrant_cost_per_unit)

        if system == "Pack"
            print("\nbattery_pack_price_to_OEM_without_BMS_Integration = ",battery_pack_price_to_OEM_without_BMS_Integration ," ")
            print("\ntotal_battery_cost_to_OEM = ",total_battery_cost_to_OEM ," ")
        else
            print("\ntotal_cell_cost_to_OEM             = ", cell_price_to_OEM)
        end

        println("\n\n************************************************ Results ************************************************")
        # println("\nBattery ", i)
        print("\nMFG Capacity                       = ", mfg_capacity, " MWh")
        print("\nCell kWh                           = ", cell_energy_storage)
        print("\nJulia dollars/kWh                  = ", dollars_per_kWh)
    end

    if breakdown
        return total_unit_cost_all_materials  ,
        total_direct_labor_per_unit           ,
        variable_overhead_per_unit            ,
        GSA_per_unit                          ,
        R_and_D_per_unit                      ,
        depreciation_per_unit                 ,
        profits_after_taxes_per_unit          ,
        warrant_cost_per_unit                 ,
        cell_energy_storage

    elseif materials_breakdown
        cathode_cost =cost_per_cell_positive_active_material+
            cost_per_cell_positive_current_collector+
            cost_per_cell_positive_conductive+
            cost_per_cell_binder_solvent_NMP+
            cost_per_cell_positive_binder

        anode_cost = cost_per_cell_negative_active_material+
            cost_per_cell_negative_binder+
            cost_per_cell_negative_conductive+
            cost_per_cell_negative_current_collector

        electrolyte_cost = cost_per_cell_electrolyte

        others_cost = cost_per_cell_cell_container+
            cost_per_cell_negative_terminal_assembly+
            cost_per_cell_positive_terminal_assembly+
            cost_per_cell_separator

        return cathode_cost,anode_cost,electrolyte_cost,others_cost
    else
        return dollars_per_kWh, mfg_capacity
    end
end

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
@testset "Utils" begin

    @testset "Test soft_theta" begin
        x = -1:0.1:1
        @test eltype(MicroscopyTools.soft_theta.(Float32.(x))) == Float32
        @test MicroscopyTools.soft_theta.(x, Ref(13.1)) ≈ 1 ./ (1 .+exp.(-2 .*13.1 .*x))
        @test MicroscopyTools.soft_theta_pw(0.01) == 1.0
        @test typeof(MicroscopyTools.soft_theta_pw(0.01f0)) == Float32
        @test MicroscopyTools.soft_theta_pw(-0.01) == 0.0
        @test MicroscopyTools.soft_theta_pw(-0.0) ≈ 0.5
        @test MicroscopyTools.soft_theta_pw(1/3.0, 1.0) ≈ 0.75
    end

    @testset "Test soft_delta" begin

        T = Float64
        f(x, k) = one(T) / (k * √(T(π))) * exp(-abs2(x / k))

        x = -1:0.1:1
        @test eltype(MicroscopyTools.soft_delta.(Float32.(x))) == Float32
        @test MicroscopyTools.soft_delta.(x, Ref(13.1)) ≈ f.(x, 13.1)

        @test typeof(MicroscopyTools.soft_delta_pw(0.01f0)) == Float32
        @test MicroscopyTools.soft_delta_pw(0.01f0) == 0.0
        @test MicroscopyTools.soft_delta_pw(-0.01) == 0.0
        @test MicroscopyTools.soft_delta_pw(-0.0) ≈ 1.0
        @test MicroscopyTools.soft_delta_pw(1/3.0, 1.0) ≈ 0.75
    end

end

ext: jl

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges

struct VariableDummyBridge <: MOIB.Variable.AbstractBridge
    id::Int
end
MOIB.Variable.unbridged_map(::VariableDummyBridge, ::MOI.VariableIndex) = nothing

map = MOIB.Variable.Map()
@testset "Empty" begin
    @test isempty(map)
    @test !MOIB.Variable.has_bridges(map)
    @test isempty(keys(map))
end

b1 = VariableDummyBridge(1)
set1 = MOI.EqualTo(0.0)
F1 = MOI.SingleVariable
S1 = typeof(set1)
v1, c1 = MOIB.Variable.add_key_for_bridge(map, () -> b1, set1)
cannot_unbridge_err = ErrorException(
    "Cannot unbridge function because some variables are bridged by variable" *
    " bridges that do not support reverse mapping, e.g., `ZerosBridge`."
)
@testset "Scalar set" begin
    @test v1.value == c1.value == -1
    @test MOIB.Variable.constraint(map, v1) == c1
    @test haskey(map, v1)
    @test map[v1] == b1
    @test MOIB.Variable.constrained_set(map, v1) == S1
    @test MOIB.Variable.number_with_set(map, S1) == 1
    @test MOIB.Variable.constraints_with_set(map, S1) == [c1]
    @test MOIB.Variable.function_for(map, c1) == MOI.SingleVariable(v1)
    @test_throws cannot_unbridge_err MOIB.Variable.unbridged_function(
        map, MOI.VariableIndex(1))

    @test MOIB.Variable.number_of_variables(map) == 1
    @test MOIB.Variable.list_of_constraint_types(map) == Set([(F1, S1)])
    @test length(map) == 1
    @test !isempty(map)
    @test MOIB.Variable.has_bridges(map)
    @test collect(keys(map)) == [v1]
    @test collect(values(map)) == [b1]
end

b2 = VariableDummyBridge(2)
set2 = MOI.Zeros(4)
F2 = MOI.VectorOfVariables
S2 = typeof(set2)
v2, c2 = MOIB.Variable.add_keys_for_bridge(map, () -> b2, set2)
@testset "Vector set" begin
    @test v2[1].value == c2.value == -2
    @test MOIB.Variable.has_keys(map, v2)
    @test !MOIB.Variable.has_keys(map, v2[4:-1:1])
    for i in 1:4
        @test MOIB.Variable.constraint(map, v2[i]) == c2
        @test haskey(map, v2[i])
        @test map[v2[i]] == b2
        @test MOIB.Variable.constrained_set(map, v2[i]) == S2
        @test MOIB.Variable.length_of_vector_of_variables(map, v2[i]) == 4
        @test MOIB.Variable.index_in_vector_of_variables(map, v2[i]) == MOIB.Variable.IndexInVector(i)
    end
    @test MOIB.Variable.number_with_set(map, S2) == 1
    @test MOIB.Variable.constraints_with_set(map, S2) == [c2]
    @test MOIB.Variable.function_for(map, c2) == MOI.VectorOfVariables(v2)
    @test_throws cannot_unbridge_err MOIB.Variable.unbridged_function(
        map, MOI.VariableIndex(1))

    @test MOIB.Variable.number_of_variables(map) == 5
    @test MOIB.Variable.list_of_constraint_types(map) == Set([(F1, S1), (F2, S2)])
    @test length(map) == 2
    @test sprint(MOIB.print_num_bridges, map) == "\nwith 2 variable bridges"
    @test !isempty(map)
    @test MOIB.Variable.has_bridges(map)
    @test collect(keys(map)) == [v1; v2]
    @test collect(values(map)) == [b1, b2]
end

b3 = VariableDummyBridge(3)
set3 = MOI.Zeros(0)
v3, c3 = MOIB.Variable.add_keys_for_bridge(map, () -> b3, set3)
@testset "Vector set of length 0" begin
    @test isempty(v3)
    @test c3.value == 0

    bridges = collect(values(map))
    @test sort([b.id for b in bridges]) == 1:2
    elements = sort(collect(map), by = el -> el.second.id)
    @test elements[1].first == v1
    @test elements[1].second == b1
    @test elements[2].first == v2[1]
    @test elements[2].second == b2
end

@testset "Delete" begin
    delete!(map, v1)
    @test MOIB.Variable.number_of_variables(map) == 4
    @test length(map) == 1
    @test sprint(MOIB.print_num_bridges, map) == "\nwith 1 variable bridge"
    @test !isempty(map)
    @test MOIB.Variable.has_bridges(map)
    elements = collect(map)
    @test elements[1].first == v2[1]
    @test elements[1].second == b2
    @test MOIB.Variable.number_with_set(map, S1) == 0
    @test MOIB.Variable.number_with_set(map, S2) == 1
    @test isempty(MOIB.Variable.constraints_with_set(map, S1))
    @test MOIB.Variable.constraints_with_set(map, S2) == [c2]
    @test MOIB.Variable.list_of_constraint_types(map) == Set([(F2, S2)])
    @test collect(keys(map)) == v2
    @test collect(values(map)) == [b2]
    @test !haskey(map, v1)
    @test MOIB.Variable.has_keys(map, v2)

    rev_v2 = reverse(v2)
    err = ArgumentError(
        "`$rev_v2` is not a valid key vector as returned by `add_keys_for_bridge`."
    )
    @test_throws err delete!(map, rev_v2)

    delete!(map, v2[3])
    left = [1, 2, 4]
    @test MOIB.Variable.number_of_variables(map) == 3
    @test length(map) == 1
    @test !isempty(map)
    @test MOIB.Variable.has_bridges(map)
    elements = collect(map)
    @test elements[1].first == v2[1]
    @test elements[1].second == b2
    @test MOIB.Variable.number_with_set(map, S1) == 0
    @test MOIB.Variable.number_with_set(map, S2) == 1
    @test isempty(MOIB.Variable.constraints_with_set(map, S1))
    @test MOIB.Variable.constraints_with_set(map, S2) == [c2]
    @test MOIB.Variable.list_of_constraint_types(map) == Set([(F2, S2)])
    @test collect(keys(map)) == v2[left]
    @test collect(values(map)) == [b2]
    @test !haskey(map, v1)
    @test !haskey(map, v2[3])
    @test !MOIB.Variable.has_keys(map, v2)
    @test MOIB.Variable.has_keys(map, v2[left])
    for (j, i) in enumerate(left)
        @test MOIB.Variable.constraint(map, v2[i]) == c2
        @test haskey(map, v2[i])
        @test map[v2[i]] == b2
        @test MOIB.Variable.constrained_set(map, v2[i]) == S2
        @test MOIB.Variable.length_of_vector_of_variables(map, v2[i]) == 3
        @test MOIB.Variable.index_in_vector_of_variables(map, v2[i]) == MOIB.Variable.IndexInVector(j)
    end
    @test MOIB.Variable.function_for(map, c2) == MOI.VectorOfVariables(v2[left])

    delete!(map, v2[1])
    left = [2, 4]
    @test MOIB.Variable.number_of_variables(map) == 2
    @test length(map) == 1
    @test !isempty(map)
    @test MOIB.Variable.has_bridges(map)
    elements = collect(map)
    @test elements[1].first == v2[1]
    @test elements[1].second == b2
    @test MOIB.Variable.number_with_set(map, S1) == 0
    @test MOIB.Variable.number_with_set(map, S2) == 1
    @test isempty(MOIB.Variable.constraints_with_set(map, S1))
    @test MOIB.Variable.constraints_with_set(map, S2) == [c2]
    @test MOIB.Variable.list_of_constraint_types(map) == Set([(F2, S2)])
    @test collect(keys(map)) == v2[left]
    @test collect(values(map)) == [b2]
    @test !haskey(map, v1)
    @test !haskey(map, v2[1])
    @test !haskey(map, v2[3])
    @test !MOIB.Variable.has_keys(map, v2)
    @test !MOIB.Variable.has_keys(map, v2[[2, 3, 4]])
    @test !MOIB.Variable.has_keys(map, v2[[1, 2, 4]])
    @test MOIB.Variable.has_keys(map, v2[left])
    for (j, i) in enumerate(left)
        @test MOIB.Variable.constraint(map, v2[i]) == c2
        @test haskey(map, v2[i])
        @test map[v2[i]] == b2
        @test MOIB.Variable.constrained_set(map, v2[i]) == S2
        @test MOIB.Variable.length_of_vector_of_variables(map, v2[i]) == 2
        @test MOIB.Variable.index_in_vector_of_variables(map, v2[i]) == MOIB.Variable.IndexInVector(j)
    end
    @test MOIB.Variable.function_for(map, c2) == MOI.VectorOfVariables(v2[left])

    delete!(map, v2[left])
    @test MOIB.Variable.number_of_variables(map) == 0
    @test length(map) == 0
    @test isempty(map)
    @test MOIB.Variable.has_bridges(map)
    @test isempty(collect(map))
    @test MOIB.Variable.number_with_set(map, S1) == 0
    @test MOIB.Variable.number_with_set(map, S2) == 0
    @test isempty(MOIB.Variable.constraints_with_set(map, S1))
    @test isempty(MOIB.Variable.constraints_with_set(map, S2))
    @test isempty(MOIB.Variable.list_of_constraint_types(map))
    @test isempty(keys(map))
    @test isempty(values(map))
    @test !haskey(map, v1)
    @test !MOIB.Variable.has_keys(map, v2)
    @test !MOIB.Variable.has_keys(map, v2[left])
end

@testset "EmptyMap" begin
    map = MOIB.Variable.EmptyMap()
    empty!(map)
    @test isempty(map)
    @test length(map) == 0
    @test isempty(keys(map))
    @test isempty(values(map))
    @test !MOIB.Variable.has_bridges(map)
    @test iszero(MOIB.Variable.number_of_variables(map))
    @test iszero(MOIB.Variable.number_with_set(map, S1))
    @test iszero(MOIB.Variable.number_with_set(map, S2))
    @test isempty(MOIB.Variable.constraints_with_set(map, S1))
    @test isempty(MOIB.Variable.constraints_with_set(map, S2))
    @test sprint(MOIB.print_num_bridges, map) == ""
end


text:
const VALID_BAR_SIZES = (
    Second(1), Second(5), Second(10), Second(15), Second(30),
    Minute(1), Minute(2), Minute(3), Minute(5), Minute(10), Minute(15), Minute(20), Minute(30),
    Hour(1), Hour(2), Hour(3), Hour(4), Hour(8),
    Day(1),
    Week(1),
    Month(1),
)

const VALID_DURATION_UNITS = (Second, Day, Week, Month, Year)

const VALID_HISTORICAL_DATA_TYPE = (
    "TRADES",
    "MIDPOINT",
    "BID",
    "ASK",
    "BID_ASK",
    "ADJUSTED_LAST",
    "HISTORICAL_VOLATILITY",
    "OPTION_IMPLIED_VOLATILITY",
    "REBATE_RATE",
    "FEE_RATE",
    "YIELD_BID",
    "YIELD_ASK",
    "YIELD_BID_ASK",
    "YIELD_LAST",
)

# TODO Timezone
function historical_data(
    ib::IB,
    asset::Asset;
    end_datetime::DateTime = now(),
    duration::TimePeriod = Week(1),
    barsize::TimePeriod = Minute(1),
    use_regular_trading_hours::Bool=true,
    datatype::String = "MIDPOINT",
)
    check_end_datetime(end_datetime)
    check_duration(dt)
    check_barsize(sz)
    check_historical_data_type(datatype)
    datas = []
    while true
        data = ib.reqHistoricalData(
            asset,
            endDateTime = end_datetime,
            durationStr = duration,
            barsizeSetting = barsize,
            whatToShow = datatype,
            useRTH = use_regular_trading_hours
        )
        @info data
        data === nothing && break
        push!(datas, data)
        @info data[1].data
        break
    end
end

function duration_str(dt::TimePeriod)
    check_duration(dt)
    unit = dt isa Second ? "S" :
        dt isa Day ? "D" :
        dt isa Week ? "W" :
        dt isa Month ? "M" :
        dt isa Year ? "Y" :
        error("This shouldn't happen")
    return "$(dt.value) $unit"
end

function barsize_str(sz::TimePeriod)
    check_barsize(sz)
    unit = sz isa Second ? "secs" :
        sz isa Minute ? (sz == Minute(1) ? "min" : "mins") :
        sz isa Hour ? (sz == Hour(1) ? "hour" : "hours") :
        sz isa Day ? "day" :
        sz isa Week ? "week" :
        sz isa Month ? "month" :
        error("This shouldn't happen")
    return "$(dt.value) $unit"
end

function check_end_datetime(t::DateTime)
    t > now() && error("End date time cannot be in the future: $t")
    return nothing
end

function check_duration(dt::TimePeriod)
    idx = findfirst(T -> typeof(dt) === T, VALID_DURATION_UNITS)
    idx === nothing && error("Invalid duration, expected unit of one of $(VALID_DURATION_UNITS)")
    return nothing
end

function check_barsize(sz::TimePeriod)
    Tlargest = typeof(first(canonicalize(sz).periods))
    sz = round(sz, Tlargest)
    idx = findfirst(x -> x == sz, VALID_BAR_SIZES)
    idx === nothing && error("Invalid dt $sz, expected one of $VALID_BAR_SIZES")
    return nothing
end

function check_historical_data_type(type::String)
    idx = findfirst(t -> t == type, VALID_HISTORICAL_DATA_TYPE)
    idx === nothing && error("Invalid historical data type, expected unit of one of $(VALID_HISTORICAL_DATA_TYPE)")
    return nothing
end


text:
using StatsBase

"""
    smd_combined = combinelocalizations(smd::SMD)

Combine clustered localizations in `smd` into higher precision localizations.
"""
function combinelocalizations(smd::SMD)
    # Isolate and sort some arrays from `smd`.
    connectID = smd.connectID
    sortindices = sortperm(connectID)
    connectID = connectID[sortindices]
    x = smd.x[sortindices]
    y = smd.y[sortindices]
    x_se = smd.x_se[sortindices]
    y_se = smd.y_se[sortindices]
    framenum = smd.framenum[sortindices]
    datasetnum = smd.datasetnum[sortindices]

    # Loop over clusters and combine their localization data as appropriate.
    nperID = counts(connectID)
    nperID = nperID[nperID .!= 0]
    ncumulative = [0; cumsum(nperID)]
    nclusters = length(nperID)
    smd_combined = deepcopy(smd)
    for nn = 1:nclusters
        indices = (1:nperID[nn]) .+ ncumulative[nn]
        smd_combined.x[nn] = StatsBase.mean(x[indices], weights(1.0/x_se[indices].^2))
        smd_combined.y[nn] = StatsBase.mean(y[indices], weights(1.0/y_se[indices].^2))
        smd_combined.x_se[nn] = sqrt(1.0 / sum(1.0./x_se[indices].^2))
        smd_combined.y_se[nn] = sqrt(1.0 / sum(1.0./y_se[indices].^2))
        smd_combined.connectID[nn] = connectID[indices[1]]
        smd_combined.framenum[nn] = framenum[indices[1]]
        smd_combined.datasetnum[nn] = datasetnum[indices[1]]
    end
    smd_combined = isolatesmd(smd_combined, 1:nclusters)

    return smd_combined
end


text:
function update_history!(s::Sampler, t::Int, theta, knots)
    s.N[t+1] = length(knots)
    s.knots[t+1,1:s.N[t+1]] = knots
    s.theta[t+1,1:(s.N[t+1]+s.prior.Q)] = theta
end

"""
    update!(s::sampler, par)

Perform one iteration of the reversibel-jump MCMC. First all the spline weights are updated and then
a birth (add new knot), death (remove existing knot) or move (move existing knot) step is attempted.
The remaining parameters `par` of the model are just passed on to the likelihood and must be updated
outside of this function.
"""
function update!(s::Sampler, par)
    t = s.t
    Nt = s.N[t]
    Kt = Nt + s.prior.Q

    ## update theta
    theta_ac = zeros(Bool, Kt)
    s.theta[t+1,:] = s.theta[t,:]

    tune_param = get_tune_param(s.tuner)
    loglik_init = NaN

    for j in 1:Kt
        s.theta[t+1,j], theta_ac[j], loglik_init = sample_theta(s, t, j, par, tune_param, loglik_init)
    end
    thac = any(theta_ac) ## 1 if any component of theta has changed, 0 otherwise
    s.accept[4] += thac

    update!(s.tuner, thac)

    ## update theta variance
    s.attempt[5] += 1
    mu = sum(s.theta[t+1,1:Kt].^2) ##transpose(theta[t+1,:])*P*theta[t+1,:]
    s.v[t+1] = rand(InverseGamma(shape(s.prior.v) + Kt/2, scale(s.prior.v) + mu/2), 1)[1]
    s.accept[5] += 1

    if s.fixed_knots
        update_history!(s, t, s.theta[t+1,1:Kt], s.knots[t,1:Nt])
    else
        ## prior probabilities of attempting birth, death and move of a knot given current number of knots
        pb = pbirth(Nt, s.prior.N_max, s.prior.N)
        pd = pdeath(Nt, s.prior.N)
        pm = 1 - pb - pd

        u = rand()

        if (u <= pb) & (Nt < s.prior.N_max)  ## attempt birth step
            s.attempt[1] += 1
            s.accept[1] += add_knot!(s, t, par)
        elseif (u > pb) & (u <= pb+pd) & (Nt > 0) ## attempt death step
            s.attempt[2] += 1
            s.accept[2] += remove_knot!(s, t, par)
        elseif (u > pb+pd) & (Nt > 0) ## attempt move step
            s.attempt[3] += 1
            s.accept[3] += move_knot!(s, t, par)
        end
    end

    s.t += 1
end

"""
    sample_theta(s::Sampler, t::Int, j::Int, par, gamma::Float64, loglik_init=NaN)

Draw a single sample from the full conditional of `theta_j` given all other parameters `par`.
The proposal is generated from a normal distribution with mean equal to `theta[t,j]` and standard
error equal to `gamma`.
"""
function sample_theta(s::Sampler, t::Int, j::Int, par, gamma::Float64, loglik_init=NaN)
    theta = s.theta[t+1, 1:(s.N[t]+s.prior.Q)]
    knots = s.prior.cand_knots[s.knots[t, 1:s.N[t]]]
    outer = s.prior.outer
    v = s.v[t]

    m = length(theta)
    d = Normal(0, gamma)
    tmp = copy(theta)

    function log_target(x::Float64)
        tmp[j] = x
        sp = spline(outer, knots, exp.(tmp))
        s.loglik(sp, par) - s.prior.theta_penalty(tmp, v)
    end

    if isnan(loglik_init)
        loglik_init = log_target(theta[j])
    end

    prop = rand(d) + theta[j]
    loglik_prop = log_target(prop)

    if log(rand()) < (loglik_prop - loglik_init)
        return prop, true, loglik_prop
    else
        return theta[j], false, loglik_init
    end
end

function pbirth(N::Int, N_max::Int, prior_N::DiscreteUnivariateDistribution)
    if N < N_max
        0.4*min(1, pdf(prior_N, N+1)/pdf(prior_N, N))
    else
        0.0
    end
end

function pdeath(N::Int, prior_N::DiscreteUnivariateDistribution)
    if N > 1
        0.4*min(1, pdf(prior_N, N-1)/pdf(prior_N, N))
    else
        0.0
    end
end

function find_interval(knots::Vector{T}, x::T) where T <: Real
    n::Int = 0
    for kn in knots
        if x < kn
            break
        end
        n += 1
    end
    n
end

"""
     add_knot(s::Sampler, t::Int, par)

Attempt to add knot randomly selected from all unoccupied candidate knots to the spline.
The spline weights stored in `s.theta[t+1,:]` are adjusted when move is accepted.
The locations of all other knots are not changed.
"""
function add_knot!(s::Sampler, t::Int, par)
    cand_knots = s.prior.cand_knots
    outer = s.prior.outer
    Q = s.prior.Q

    ## randomly select unoccupied candidate knot
    cindex = sample((1:length(cand_knots))[.!s.knot_ocp])

    knots = s.knots[t,1:s.N[t]]
    theta = s.theta[t+1,1:(s.N[t]+s.prior.Q)]
    v = s.v[t+1]

    ckci = cand_knots[cindex]
    ckk = cand_knots[knots]

    ## find index of left endpoint of interval containing candidate knot
    n = find_interval(ckk, ckci)
    j = n+Q

    ## new knots
    new_xi = copy(knots)
    insert!(new_xi, n+1, cindex)

    ## new control points
    new_theta = copy(theta)
    insert!(new_theta, n+2, 0.0)

    pk = Spline.pad_knots(ckk, outer, Q-1)

    ## @assert Q > 2
    a = (j-Q+2):(j-1) ##(n+2):(n+Q-1)
    r = (ckci .- pk[a]) ./ (pk[a.+(Q-1)] - pk[a])
    new_theta[a] = log.(r.*exp.(theta[a]) + (1 .- r).*exp.(theta[a.-1]))

    u = rand(1)[1]
    ##u = (cand_knots[cindex] - pk[j]) ./ (pk[j+(Q-1)] - pk[j])
    new_theta[j] = log(u*exp(theta[j]) + (1-u)*exp(theta[j-1]))

    ## acceptance probability
    logRPT = -log(2*pi*v)/2 + (sum(theta.^2) - sum(new_theta.^2))/(2*v)

    sp = spline(outer, ckk, exp.(theta))
    new_sp = spline(outer, cand_knots[new_xi], exp.(new_theta))
    logRL = s.loglik(new_sp, par) - s.loglik(sp, par)

    ## Jacobian
    J = (exp(theta[j]) - exp(theta[j-1])) / exp(new_theta[j])
    for i in (j-Q+2):(j-1)
        r = (ckci - pk[i]) / (pk[i+Q-1] - pk[i])
        J *= r*exp(theta[i] - new_theta[i])
    end

    if log(rand(1)[1]) < (logRPT + logRL + log(abs(J)))
        update_history!(s, t, new_theta, new_xi)
        s.knot_ocp[cindex] = true
        true
    else
        update_history!(s, t, theta, knots)
        false
    end
end

"""
    remove_knot(s::Sampler, t::Int, n::Int, par)

Attempt to remove knot at location `n` from the spline.
The spline weights stored in `s.theta[t+1,:]` are adjusted when move is accepted.
The locations of all other knots are not changed.
"""
function remove_knot!(s::Sampler, t::Int, par)
    Q = s.prior.Q
    outer = s.prior.outer
    cand_knots = s.prior.cand_knots

    v = s.v[t+1]
    theta = s.theta[t+1, 1:(s.N[t]+s.prior.Q)]
    knots = s.knots[t, 1:s.N[t]]

    n = sample(1:s.N[t])
    cindex = knots[n]

    ## new knots
    new_xi = copy(knots)
    deleteat!(new_xi, n)

    j = n + Q

    ## new control points
    new_theta = copy(theta)
    deleteat!(new_theta, j-1)

    pk = Spline.pad_knots(cand_knots[new_xi], outer, Q-1)

    ## @assert Q > 2
    for i in (j-Q+1):(j-2)
        r = (cand_knots[cindex] .- pk[i]) / (pk[i+Q-1] - pk[i])
        nt_prop = exp(new_theta[i])/r .- exp(new_theta[i-1])*(1 - r)/r

        if nt_prop > 0
            new_theta[i] = log(nt_prop)
        else
            new_theta[i] = -100
        end
    end

    ## Jacobian
    J = exp(new_theta[j-1])/(exp(theta[j-1]) - exp(theta[j-2]))
    for i in (j - Q + 2):(j-2)
        r = (cand_knots[cindex] - pk[i]) / (pk[i+Q-1] - pk[i])
        J *= exp(new_theta[i]) / (r*exp(theta[i]))
    end

    ## acceptance probability
    logRPT = log(2*pi*v)/2 + (sum(theta.^2) - sum(new_theta.^2))/(2*v)

    sp = spline(outer, cand_knots[knots], exp.(theta))
    new_sp = spline(outer, cand_knots[new_xi], exp.(new_theta))
    logRL = s.loglik(new_sp, par) - s.loglik(sp, par)

    if log(rand(1)[1]) < (logRPT + logRL + log(abs(J)))
        update_history!(s, t, new_theta, new_xi)
        s.knot_ocp[cindex] = false
        true
    else
        update_history!(s, t, theta, knots)
        false
    end
end

"""
    move_knot(s::Sampler, t::Int, par)

Attempt to move randomly selected knot to a new location randomly chosen from all
unoccupied neighbouring knots (all knots between the selected knot and the next occupied knots to the left and right).
The ordering of the occupied knots is not changed by the move.
The spline weights stored in `s.theta[t+1,:]` are adjusted when move is accepted.
"""
function move_knot!(s::Sampler, t::Int, par)
    cand_knots = s.prior.cand_knots
    outer = s.prior.outer

    Nt = s.N[t]
    Kt = Nt + s.prior.Q

    s.N[t+1] = Nt

    knots = s.knots[t,1:Nt]
    theta = s.theta[t+1,1:Kt]

    n = sample(1:Nt)
    cindex = knots[n]
    s.knots[t+1,1:Nt] = s.knots[t,1:Nt]

    ## get range of unoccupied neighbouring knots (including the currently selected knot to be moved)
    if (n == 1) & (n < Nt)
        nbg_knots = 1:(knots[n+1]-1)
    elseif (n == length(knots)) & (n > 1)
        nbg_knots = (knots[n-1]+1):length(cand_knots)
    elseif (n > 1) & (n < Nt)
        nbg_knots = (knots[n-1]+1):(knots[n+1]-1)
    else
        nbg_knots = 1:length(cand_knots)
    end

    ## check if all neighbouring candidate knots already occupied
    if length(nbg_knots) <= 1
        return false
    end

    ## uniformly sample new position from all candidate knots between bounds
    new_cindex = sample(nbg_knots)

    if new_cindex == cindex ## new position same as old
        return false
    else
        tmp = cand_knots[knots]
        tmp[n] = cand_knots[new_cindex]

        sp = spline(outer, cand_knots[knots], exp.(theta))
        new_sp = spline(outer, tmp, exp.(theta))
        logRL = s.loglik(new_sp, par) - s.loglik(sp, par)

        if log(rand(1)[1]) < logRL
            ## after the moved knot is still at position `n` among all occupied knots
            s.knots[t+1,n] = new_cindex
            s.knot_ocp[cindex] = false
            s.knot_ocp[new_cindex] = true
            true
        else
            false
        end
    end
end


text:
"""
    host2device_bandwidth([memsize::UnitPrefixedBytes=GiB(0.5)]; kwargs...)

Performs a host-to-device memory copy benchmark (time measurement) and
returns the host-to-device bandwidth estimate (in GiB/s) derived from it.

**Keyword arguments:**
* `nbench` (default: `10`): number of time measurements (i.e. p2p memcopies)
* `verbose` (default: `true`): set to false to turn off any printing.
* `stats` (default: `false`): when `true` shows statistical information about the benchmark.
* `times` (default: `false`): toggle printing of measured times.
* `dtype` (default: `Cchar`): used data type.
* `io` (default: `stdout`): set the stream where the results should be printed.

**Examples:**
```julia
host2device_bandwidth()
host2device_bandwidth(MiB(1024))
host2device_bandwidth(KiB(20_000); dtype=Int32)
```
"""
function host2device_bandwidth(
    memsize::UnitPrefixedBytes=GiB(0.5);
    dtype=Cchar,
    DtoDfactor=true,
    verbose=true,
    io::IO=stdout,
    kwargs...,
)
    N = Int(bytes(memsize) ÷ sizeof(dtype))
    mem_host = rand(dtype, N)
    mem_host_pinned = Mem.pin(rand(dtype, N))
    mem_gpu = CUDA.rand(dtype, N)
    # mem_gpu2 = CUDA.rand(dtype, N)

    # if verbose
    #     gpu = device(mem_gpu)
    #     println("Memsize: $(Base.format_bytes(sizeof(mem_gpu)))")
    #     println("GPU: ", gpu, " - ", name(gpu), "\n")
    # end

    _perform_memcpy(mem_host, mem_gpu; title="Host <-> Device", verbose, io=io, kwargs...)
    verbose && println(io)
    _perform_memcpy(
        mem_host_pinned, mem_gpu; title="Host (pinned) <-> Device", verbose, io=io, kwargs...
    )
    # verbose && println()
    # _perform_memcpy(mem_gpu, mem_gpu2; title="Device <-> Device (same device)", DtoDfactor, verbose, kwargs...)
    return nothing
end

function _perform_memcpy(
    mem1,
    mem2;
    title="",
    nbench=10,
    times=false,
    stats=false,
    DtoDfactor=false,
    verbose=true,
    io::IO=stdout
)
    NVTX.@range "host2dev: $title" begin
        sizeof(mem1) == sizeof(mem2) || error("sizeof(mem1) != sizeof(mem2)")
        ts = zeros(nbench)
        NVTX.@range "host2dev: bench loop" begin
            @inbounds for i in 1:nbench
                if i % 2 == 0
                    ts[i] = CUDA.@elapsed copyto!(mem1, mem2)
                else
                    ts[i] = CUDA.@elapsed copyto!(mem2, mem1)
                end
            end
        end

        t_min = minimum(ts)
        t_max = maximum(ts)
        t_avg = mean(ts)

        actual_memsize_GiB = sizeof(mem1) * 2^(-30)
        if DtoDfactor
            actual_memsize_GiB *= 2 # must count both the read and the write here (taken from p2pBandwidthLatencyTest cuda sample....)
        end
        bws = actual_memsize_GiB ./ ts
        bw_min = minimum(bws)
        bw_max = maximum(bws)
        bw_avg = mean(bws)

        if verbose
            if times
                println(io,"t_min: $t_min")
                println(io,"t_max: $t_max")
                println(io,"t_avg: $t_avg")
            end
            printstyled(io,"$(title) Bandwidth (GiB/s):\n"; bold=true)
            if stats
                print(io, " ├ max: ")
                printstyled(io, round(bw_max; digits=2), "\n"; color=:green, bold=true)
                println(io, " ├ min: ", round(bw_min; digits=2))
                println(io, " ├ avg: ", round(bw_avg; digits=2))
                print(io," └ std_dev: ")
                printstyled(io,round(std(bws); digits=2), "\n"; color=:yellow, bold=true)
            else
                print(io," └ max: ")
                printstyled(io,round(bw_max; digits=2), "\n"; color=:green, bold=true)
            end
        end
    end
    return bw_max
end


text:
import Combinatorics
import StatsModels

"""
"""
function generate_formula(lhs::AbstractVector{<:StatsModels.AbstractTerm},
                          rhs::AbstractVector{<:StatsModels.AbstractTerm})
    lhs_deepcopy = deepcopy(lhs)
    rhs_deepcopy = deepcopy(rhs)
    lhs_sum = sum(lhs_deepcopy)
    rhs_sum = sum(rhs_deepcopy)
    result = lhs_sum ~ rhs_sum
    return result
end

"""
"""
function generate_formula(
        dependent_variable::Symbol,
        independent_variables::AbstractVector{<:Symbol};
        intercept::Bool = true,
        interactions::Integer = 1,
        )
    result = generate_formula(
        [dependent_variable],
        independent_variables;
        intercept = intercept,
        interactions = interactions,
        )
    return result
end

"""
"""
function generate_formula(
        dependent_variables::AbstractVector{<:Symbol},
        independent_variables::AbstractVector{<:Symbol};
        intercept::Bool = true,
        interactions::Integer = 1,
        )
    if length(dependent_variables) < 1
        error("length(dependent_variables) must be >= 1")
    end
    if length(independent_variables) < 1
        error("length(independent_variables) must be >= 1")
    end
    if interactions < 1
        error("interactions must be >= 1")
    end
    if interactions > length(independent_variables)
        error("interactions must be <= the number of independent variable")
    end
    if intercept
        intercept_term = StatsModels.term(1)
    else
        intercept_term = StatsModels.term(0)
    end
    lhs_terms = Vector{StatsModels.AbstractTerm}(undef, 0)
    for dep_var in deepcopy(dependent_variables)
        push!(lhs_terms, StatsModels.term(dep_var))
    end
    rhs_terms = Vector{StatsModels.AbstractTerm}(undef, 0)
    push!(rhs_terms, intercept_term)
    for interaction_term in generate_interaction_terms_up_to_level(
            deepcopy(independent_variables),
            interactions)
        push!(rhs_terms, interaction_term)
    end
    result = generate_formula(lhs_terms, rhs_terms)
    return result
end

function generate_interaction_terms_up_to_level(
        xs::AbstractVector{<:Symbol},
        interactions::Integer,
        )
    all_interaction_terms = Vector{StatsModels.AbstractTerm}(undef, 0)
    for level = 1:interactions
        append!(all_interaction_terms,
                generate_interaction_terms_at_level(xs,level))
    end
    unique!(all_interaction_terms)
    return all_interaction_terms
end

function generate_interaction_terms_at_level(
        xs::AbstractVector{<:Symbol},
        level::Integer,
        )
    interaction_sets = collect(
        Combinatorics.combinations(
            xs,
            level,
            )
        )
    interaction_terms = Vector{StatsModels.AbstractTerm}(undef, 0)
    for interaction_set in interaction_sets
        push!(interaction_terms,
              reduce(&, StatsModels.term.(interaction_set)))
    end
    unique!(interaction_terms)
    return interaction_terms
end


text:
using RayTracer, Zygote, Flux, Images, Statistics

mkdir("images")

function create_and_save(color, val)
    img = get_image(color, screen_size.w, screen_size.h)
    save("images/triangle_$(val).jpg", img)
end

screen_size = (w = 200, h = 200)

light = PointLight(Vec3(1.0), 20000.0, Vec3(5.0, 5.0, -10.0))

light_perturbed = PointLight(Vec3(1.0), 20000.0, Vec3(1.0, 2.0, -7.0))

eye_pos = Vec3(0.0, 0.0, -5.0)

scene = [
    Triangle(Vec3(-1.7, 1.0, 0.0), Vec3(1.0, 1.0, 0.0), Vec3(-0.5, -1.0, 0.0),
             color = rgb(0.0, 1.0, 0.0), reflection = 0.5)
    ]

origin, direction = get_primary_rays(Float64, screen_size.w, screen_size.h, 45, eye_pos);

color = raytrace(origin, direction, scene, light, eye_pos, 0)

create_and_save(color, "original")

color_guess = raytrace(origin, direction, scene, light_perturbed, eye_pos, 0)

create_and_save(color_guess, "initial_guess")

function process_image(im::Vec3{T}, width, height) where {T}
    im_arr = reshape(hcat(im.x, im.y, im.z), (width, height, 3))

    return reshape(im_arr, (width, height, 3, 1))
end

function loss_fn(θ, img)
    rendered_color = raytrace(origin, direction, scene, θ, eye_pos, 0)
    rendered_img = process_image(rendered_color, screen_size.w, screen_size.h)
    loss = mean(abs2.(rendered_img .- img))
    @show loss
    return loss
end

image = process_image(color, screen_size.w, screen_size.h);

opt = ADAM(0.05)

for i in 1:9000
    gs = gradient(x -> loss_fn(x, image), light_perturbed)[1]
    update!(opt, light_perturbed.position, gs.position)
    if i % 5 == 0
        @info "$i iterations completed"
        create_and_save(raytrace(origin, direction, scene, light_perturbed, eye_pos, 0), i)
        display(light_perturbed)
    end
end

@info "Light Position Optimized"

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
using ComplexSystems
using Plots
gr()

fl65 = Dict{String,String}("F" => "F[-F]F[+F]F")
HW = mult_evol_lsystem("F",fl65,3)
splitted = split(HW,"");
plot();
let
    global x = 0.
    global y = 0.
    global a = 0.
    global p
    global variables=[(x,y,a)]
    Δ = 25
    global da = Δ/180*π
    for s in splitted
        if s == "F" || s == "G"
            p,x,y=turtle_step(x,y,s,a)
        elseif s == "+" || s == "-"
            a = update_angle(a,da,1,s)
        elseif s == "["
            push!(variables,(x,y,a))
            x,y,a = variables[end]
        elseif s == "]"
            x,y,a = variables[end]
            pop!(variables)
        end
        p
    end
end

mf = Dict{String,String}("F" => "FF[+F][-F]","+" =>
"+F[+2F][-2F]","-"=>"-F[-3F][+3F]")
HW2 = mult_evol_lsystem("F",mf,4);
splitted = split(HW2,"");

plot();
let
    global x = 0.
    global y = 0.
    global a = 0.
    global p
    global variables=[(x,y,a)]
    Δ = 45
    global da = Δ/180*π
    for s in splitted
        if s == "F" || s == "G"
            p,x,y=turtle_step(x,y,s,a)
        elseif s == "+" || s == "-"
            a = update_angle(a,da,1,s)
        elseif s == "["
            push!(variables,(x,y,a))
            x,y,a = variables[end]
        elseif s == "]"
            x,y,a = variables[end]
            pop!(variables)
        end
        p
    end
end


text:
export Beam

mutable struct Beam <: AbstractElement
    id::String
    hid::Int
    node1::Node
    node2::Node
    material::Material
    section::BeamSection

    release::Vector{Bool}

    elm_type::String
    mass_type::String

    center::Vector{Float64}
    l::Float64
    T::Matrix{Float64} #transform_matrix
end

function Beam(id,hid,node1,node2,material,section;elm_type="eular_shear",mass_type="concentrate")
    tol=1e-6
    o=node1.loc
    pt1=node2.loc
    pt2=node1.loc
    if abs(pt2[1]-pt1[1])<tol && abs(pt2[2]-pt1[2])<tol
        pt2=pt2.+[1,0,0]
    else
        pt2=pt2.+[0,0,1]
    end
    csys=CSys(o,pt1,pt2)
    T=zeros(12,12)
    T[1:3,1:3]=T[4:6,4:6]=T[7:9,7:9]=T[10:12,10:12]=csys.T
    l=norm(node1.loc-node2.loc)
    release=Bool.(zeros(12))
    Beam(string(id),hid,node1,node2,material,section,release,elm_type,mass_type,(pt1+pt2)/2,l,T)
end

for (root,dirs,files) in walkdir(joinpath(@__DIR__,"beams"))
    for file in files
        if file[end-2:end]==".jl"
            include(joinpath(@__DIR__,"beams",file))
        end
    end
end

function integrateK(beam::Beam)::Matrix{Float64}
    if beam.elm_type=="eular_shear"
        return K_eular_shear(beam::Beam)
    end
end

function integrateKσ(beam::Beam,σ::Vector{Float64})::Matrix{Float64}
    if beam.elm_type=="eular_shear"
        return K2_eular_shear(beam::Beam)
    end
end

function static_condensation(K,P,rDOF::Vector{Int})
    if isempty(rDOF)
        return K,P
    end
    i=[!(x in rDOF) for x in 1:12]
    j=[(x in rDOF) for x in 1:12]
    Kᵢᵢ=K[i,i]
    Kᵢⱼ=K[i,j]
    Kⱼᵢ=K[j,i]
    Kⱼⱼ=K[j,j]
    Kⱼⱼ⁻¹=inv(Kⱼⱼ)
    Pⱼ=P[j]
    Pᵢ=P[i]
    Kᶜ=zero(K)
    Pᶜ=zero(P)
    Kᶜ[i,i]=Kᵢᵢ-Kᵢⱼ*Kⱼⱼ⁻¹*Kⱼᵢ
    Pᶜ[i]=Pᵢ-Kᵢⱼ*Kⱼⱼ⁻¹*Pⱼ
    return Kᶜ,Pᶜ
end

function integrateM(beam::Beam)::Matrix{Float64}
    E,ν=beam.material.E,beam.material.ν
    A,I₂,I₃,J,l=beam.section.A,beam.section.I₂,beam.section.I₃,beam.section.J,beam.l
    ρ=beam.material.ρ
    if beam.mass_type=="concentrate"
        return Matrix(I,12,12)*12*ρ*A*l/2
    elseif beam.mass_type=="coordinate"
        return Matrix(I,12,12)*12*ρ*A*l/2
    end
end

function calc_force(beam,beamforce)
    f₁,f₂=beamforce.f[1:6],beamforce.f[7:12]
    l=beam.l
    fi,v2i,v3i,ti,m2i,m3i,fj,v2j,v3j,tj,m2j,m3j=beamforce.f
    Nᵀf(x)=[
          0;
    (1 - 3*x^2 + 2*x^3)*(v2i + x*(-v2i + v2j));
    (1 - 3*x^2 + 2*x^3)*(v3i + x*(-v3i + v3j));
          0;
    l*(x - 2*x^2 + x^3)*(v3i + x*(-v3i + v3j));
    l*(x - 2*x^2 + x^3)*(v2i + x*(-v2i + v2j));
          0;
        (3*x^2 - 2*x^3)*(v2i + x*(-v2i + v2j));
        (3*x^2 - 2*x^3)*(v3i + x*(-v3i + v3j));
         0;
         l*(-x^2 + x^3)*(v3i + x*(-v3i + v3j));
         l*(-x^2 + x^3)*(v2i + x*(-v2i + v2j));
   ]*l

   Nᵀf2(x)=[
         -0.5*(-1 + x)*(fi + (1+x)/2*(-fi + fj));0;0;
         -0.5*(-1 + x)*(ti + (1+x)/2*(-ti + tj));0;0;
           0.5*(1 + x)*(fi + (1+x)/2*(-fi + fj));0;0;
           0.5*(1 + x)*(ti + (1+x)/2*(-ti + tj));0;0
  ]*l/2

    Pᵉ_f=hquadrature(Nᵀf,0,1)[1]+hquadrature(Nᵀf2,-1,1)[1]#Pᵉf=∫NᵀfdV #体积力

    Pᵉ_s=beamforce.s# Pᵉs=∫NᵀTdS #边界力

    f₁,f₂=beamforce.σ₀[1],beamforce.σ₀[2]
    Bᵀσ₀(x)=[-0.5*(f₁+0.5*(x+1)*(-f₁+f₂));
              0.5*(f₁+0.5*(x+1)*(-f₁+f₂))]*beam.l/2
    a,b=hquadrature(Bᵀσ₀,-1,1)[1]# Pᵉσ₀=-∫Bᵀσ₀dV #初应力
    Pᵉ_σ₀=[a,0,0,0,0,0,b,0,0,0,0,0]

    f₁,f₂=beamforce.ϵ₀[1],beamforce.ϵ₀[2]
    BᵀDϵ₀(x)=[-0.5*(f₁+0.5*(x+1)*(-f₁+f₂));
               0.5*(f₁+0.5*(x+1)*(-f₁+f₂))]*beam.l/2
    a,b=hquadrature(BᵀDϵ₀,-1,1)[1]*beam.material.E*beam.section.A# Pᵉϵ₀=∫BᵀDϵ₀dV #初应变
    Pᵉ_ϵ₀=[a,0,0,0,0,0,b,0,0,0,0,0]
    return Pᵉ_f,Pᵉ_s,Pᵉ_σ₀,Pᵉ_ϵ₀
end

function integrateP(beam::Beam,beam_force)::Vector{Float64}
    Pfᵉ,Psᵉ,Pσ₀ᵉ,Pϵ₀ᵉ=calc_force(beam,beam_force) #Pᵉ=Pᵉf+Pᵉs+Pᵉσ₀+Pᵉϵ₀
    return reshape(Pfᵉ+Psᵉ+Pσ₀ᵉ+Pϵ₀ᵉ,12)
end

# end


text:
import Base: show, showcompact

"""

# Available top level Method

### Method
```julia
*  Sample::Method             : Sampling
*  Optimize::Method           : Optimization
*  Diagnose::Method           : Diagnostics
*  Variational::Method        : Variational Bayes
```
"""
abstract type Method end

DataDict = Dict{String, Any}

mutable struct Random
  seed::Int64
end
Random(;seed::Number=-1) = Random(seed)

mutable struct Output
  file::String
  diagnostic_file::String
  refresh::Int64
end
Output(;file::String="", diagnostic_file::String="", refresh::Number=100) =
  Output(file, diagnostic_file, refresh)

mutable struct Stanmodel
  name::String
  nchains::Int
  num_warmup::Int
  num_samples::Int
  thin::Int
  id::Int
  model::String
  model_file::String
  monitors::Vector{String}
  data::Vector{DataDict}
  data_file::String
  command::Vector{Base.AbstractCmd}
  method::Method
  random::Random
  init::Vector{DataDict}
  init_file::String
  output::Output
  tmpdir::String
  useMamba::Bool
  mambaThinning::Int
end

"""
# Method Stanmodel

Create a Stanmodel.

### Constructors
```julia
Stanmodel(
  method=Sample();
  name="noname",
  nchains=4,
  num_warmup=1000,
  num_samples=1000,
  thin=1,
  model="",
  monitors=String[],
  data=DataDict[],
  random=Random(),
  init=DataDict[],
  output=Output(),
  pdir::String=pwd(),
  useMamba=true,
  mambaThinning=1
)

```
### Required arguments
```julia
* `method::Method`            : See ?Method
```

### Optional arguments
```julia
* `name::String`               : Name for the model
* `nchains::Int`               : Number of chains, if possible execute in parallel
* `num_warmup::Int`            : Number of samples used for num_warmupation
* `num_samples::Int`           : Sample iterations
* `thin::Int`                  : Stan thinning factor
* `model::String`              : Stan program source
* `data::DataDict[]`           : Observed input data as an array of Dicts
* `random::Random`             : Random seed settings
* `init::DataDict[]`           : Initial values for parameters in parameter block
* `output::Output`             : File output options
* `pdir::String`               : Working directory
* `monitors::String[] `        : Filter for variables used in Mamba post-processing
* `useMamba::Bool`             : Use Mamba Chains for diagnostics and graphics
* `mambaThinning::Int`         : Additional thinning factor in Mamba Chains
```

### Example
```julia
bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
"

stanmodel = Stanmodel(num_samples=1200, thin=2, name="bernoulli", model=bernoullimodel);
```

### Related help
```julia
?stan                          : Run a Stanmodel
?Sample                        : Sampling settings
?Method                       : List of available methods
?Output                        : Output file settings
?DataDict                      : Input data dictionaries, will be converted to R syntax
```
"""
function Stanmodel(
  method=Sample();
  name="noname",
  nchains=4,
  num_warmup=1000,
  num_samples=1000,
  thin=1,
  model="",
  monitors=String[],
  data=DataDict[],
  random=Random(),
  init=DataDict[],
  output=Output(),
  pdir::String=pwd(),
  useMamba=true,
  mambaThinning=1)

  cd(pdir)

  if useMamba
    res = ""
    try
    	res = isdefined(Main, :Mamba)
    catch
      println("Package Mamba.jl not installed, please run Pkg.add(\"Mamba\")")
      return
    end

    res = ""
    try
    	res = isdefined(Main, :Gadfly)
    catch
      println("Package Gadfly.jl not installed, please run Pkg.add(\"Gadfly\")")
      return
    end

    if typeof(res) == VersionNumber
      eval(quote
        using Mamba
        using Gadfly
      end)
    end
  end

  tmpdir = joinpath(pdir, "tmp")
  if !isdir(tmpdir)
    mkdir(tmpdir)
  end

  model_file = "$(name).stan"
  if length(model) > 0
    update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))
  end

  id::Int=0
  data_file::String=""
  init_file::String=""
  cmdarray = fill(``, nchains)

  if num_samples != 1000
    method.num_samples=num_samples
  end

  if num_warmup != 1000
    method.num_warmup=num_warmup
  end

  if thin != 1
    method.thin=thin
  end

  Stanmodel(name, nchains,
    num_warmup, num_samples, thin,
    id, model, model_file, monitors,
    data, data_file, cmdarray, method, random,
    init, init_file, output, tmpdir,
    useMamba, mambaThinning);
end

"""

# Method update_model_file

Update Stan language model file if necessary

### Method
```julia
update_model_file(
  file::String,
  str::String
)
```
### Required arguments
```julia
* `file::String`                : File holding existing Stan model
* `str::String`                 : Stan model string
```

### Related help
```julia
?Stan.Stanmodel                 : Create a StanModel
```
"""
function update_model_file(file::AbstractString, str::AbstractString)
  str2 = ""
  if isfile(file)
    resfile = open(file, "r")
    str2 = read(resfile, String)
    str != str2 && rm(file)
  end
  if str != str2
    println("\nFile $(file) will be updated.\n")
    strmout = open(file, "w")
    write(strmout, str)
    close(strmout)
  end
end

function model_show(io::IO, m::Stanmodel, compact::Bool)
  println("  name =                    \"$(m.name)\"")
  println("  nchains =                 $(m.nchains)")
  println("  num_samples =             $(m.num_samples)")
  println("  num_warmup =                   $(m.num_warmup)")
  println("  thin =                    $(m.thin)")
  println("  useMamba =                $(m.useMamba)")
  println("  mambaThinning =           $(m.mambaThinning)")
  println("  monitors =                $(m.monitors)")
  println("  model_file =              \"$(m.model_file)\"")
  println("  data_file =               \"$(m.data_file)\"")
  println("  output =                  Output()")
  println("    file =                    \"$(m.output.file)\"")
  println("    diagnostics_file =        \"$(m.output.diagnostic_file)\"")
  println("    refresh =                 $(m.output.refresh)")
  if isa(m.method, Sample)
    sample_show(io, m.method, compact)
  elseif isa(m.method, Optimize)
    optimize_show(io, m.method, compact)
  elseif isa(m.method, Variational)
    variational_show(io, m.method, compact)
  else
    diagnose_show(io, m.method, compact)
  end
end

show(io::IO, m::Stanmodel) = model_show(io, m, false)


text:
export PeronaMalik, PM1, PM2, denoise

using ImageFiltering

@qstruct PeronaMalik{PM,S}(;
    lambda::Float64=0.1,
    variant::PM = PM2(10),
    border::Symbol = :replicate,
    niter::Int=100,
    step::S=generic_step!,
    )

function generic_step!(out::AbstractArray,
                       bc::BorderArray,
                       bimg::BorderArray,
                       alg::PeronaMalik)
    cstencil = CStencil(alg.variant)
    pmstencil = PMStencil(alg.lambda)
    mapstencil!(cstencil ,  bc.inner  , bimg,   )
    mapstencil!(pmstencil, out, bimg, bc)
    out
end

function mapstencil_inner!(stencil, out, barrs...; axes)
    arrs = map(barr -> barr.inner, barrs)
    @simd for i in CartesianIndices(axes)
        @inbounds out[i] = stencil(arrs..., i)
    end
    out
end

function mapstencil_border!(stencil, out, barrs...)
    outeraxes = axes(out)
    inneraxes = map(outeraxes) do r
        r[2:end-1]
    end
    for i in EdgeIterator(outeraxes, inneraxes)
        out[i] = stencil(barrs..., i)
    end
    out
end

function mapstencil!(stencil::F, out, barrs...) where {F}
    inneraxes = map(axes(out)) do r
        r[2:end-1]
    end
    mapstencil_inner!(stencil, out, barrs..., axes=inneraxes)
    mapstencil_border!(stencil, out, barrs...)
end

function denoise(img, alg::PeronaMalik=PeronaMalik())
    if alg.variant.k isa AbstractArray
        @argcheck axes(img) == axes(alg.variant.k)
    end
    x = first(img)
    T = typeof(exp((x - x)^2 / 2))
    border = get_border(img, alg)
    out_buf = BorderArray(similar(img, T), border)
    img_buf = BorderArray(similar(img, T), border)
    c_buf   = BorderArray(similar(img, T), border)
    copy!(img_buf.inner, img)
    for i in 1:alg.niter
        alg.step(out_buf.inner, c_buf, img_buf, alg)
        out_buf, img_buf = img_buf, out_buf
    end
    img_buf.inner
end

function get_border(img::AbstractArray{T,N}, alg::PeronaMalik) where {T, N}
    dims = ntuple(_ -> 1, Val(N))
    Pad{N}(alg.border, dims, dims)
end

struct PM1{K}
    k::K
end

struct PM2{K}
    k::K
end

struct CStencil{PM}
    pm::PM
end

const PM{K} = Union{PM1{K}, PM2{K}}

@inline function get_K(pm::PM{<:Number}, index...)
    pm.k
end

@inline function get_K(pm::PM{<:AbstractArray}, index...)
    @inbounds pm.k[index...]
end

@inline function c_from_grad2(pm::PM1, K, ∇img2)
    exp(-∇img2/(K^2))
end

@inline function c_from_grad2(pm::PM2, K, ∇img2)
    1 / (1 + ∇img2/K^2)
end

@inline function (f::CStencil)(img, index)
    pm = f.pm
    K = get_K(pm, index)
    ∇img2 = dot_grad(img, img, index)
    c_from_grad2(pm, K, ∇img2)
end

function inlineit(ex)
    ret = gensym("ret")
    quote
        $(Expr(:meta, :inline))
        $(Expr(:inbounds, true))
        local $ret = $ex
        $(Expr(:inbounds, :pop))
        $ret
    end
end

"""
    dot_grad(x,y,index)

Compute (∇x⋅∇y)[index].
"""
@generated function dot_grad(x::AbstractArray{T, ndims}, y, index::CartesianIndex) where {T, ndims}
    expr_dot_grad(ndims=ndims) |> inlineit
end

function expr_dot_grad(;ndims)
    # args: x,y,index
    ret = quote
        ret = zero(Float64)
        x0 = x[index]
        y0 = y[index]
    end
    for i in 1:ndims
        index₊ = expr_shift_index(:index, ndims=ndims, i=i, shift=:(1))
        ex = quote
            index₊ = $(index₊)
            ret += (x[index₊] - x0)*(y[index₊] - y0)
        end
        append!(ret.args, ex.args)
    end
    push!(ret.args, :ret)
    ret
end

@generated function laplace(img::AbstractArray{T,ndims}, index) where {T,ndims}
    expr_laplace(ndims=ndims) |> inlineit
end

function expr_laplace(;ndims)

    ret = quote
        ret = -2*$ndims * img[index]
    end
    for i in 1:ndims
        index₊ = expr_shift_index(:index, ndims=ndims, i=i, shift=:(1))
        index₋ = expr_shift_index(:index, ndims=ndims, i=i, shift=:(-1))
        ex = quote
            index₊ = $(index₊)
            index₋ = $(index₋)
            ret += img[index₊]
            ret += img[index₋]
        end
        append!(ret.args, ex.args)
    end
    push!(ret.args, :ret)
    ret
end

function expr_shift_index(index; ndims, i, shift)
    @argcheck 1 <= i <= ndims
    args = map(1:ndims) do j
        if i == j
            :($index[$i] + $shift)
        else
            :($index[$j])
        end
    end
    :(@inbounds CartesianIndex{$ndims}($(args...)))
end

struct PMStencil
    lambda::Float64
end

@inline function (f::PMStencil)(img, c, index)
    @inbounds dudt = (dot_grad(c, img, index) + c[index] * laplace(img, index))
    @inbounds img[index] + f.lambda * dudt
end


text:
export StackFrames, ResizeImage

using ImageTransformations: imresize!

"""
    ResizeImage(img::Array{T, N})
    ResizeImage(dims::Int...) -> ResizeImage(Float32, dims...)
    ResizeImage(T::Type{<:Number}, dims::Int...)

Using BSpline method to resize the `state` field of an observation to size of `img` (or `dims`).
"""
struct ResizeImage{T,N}
    img::Array{T,N}
end

ResizeImage(dims::Int...) = ResizeImage(Float32, dims...)
ResizeImage(T::Type{<:Number}, dims::Int...) = ResizeImage(Array{T}(undef, dims))

function (p::ResizeImage)(state::AbstractArray)
    imresize!(p.img, state)
    p.img
end

"""
    StackFrames(::Type{T}=Float32, d::Int...)

Use a pre-initialized [`CircularArrayBuffer`](@ref) to store the latest several states specified by `d`. Before processing any observation, the buffer is filled with `zero{T}`.
"""
struct StackFrames{T,N}
    buffer::CircularArrayBuffer{T,N}
end

StackFrames(d::Int...) = StackFrames(Float32, d...)

function StackFrames(::Type{T}, d::Vararg{Int,N}) where {T,N}
    p = StackFrames(CircularArrayBuffer{T}(d...))
    for _ in 1:capacity(p.buffer)
        push!(p.buffer, zeros(T, size(p.buffer)[1:N-1]))
    end
    p
end

function (p::StackFrames{T,N})(state::AbstractArray) where {T,N}
    push!(p.buffer, state)
    p.buffer
end

# !!! side effect?
function Base.push!(
    cb::CircularArrayBuffer{T,N},
    stacked_data::CircularArrayBuffer{T,N},
) where {T,N}
    push!(cb, select_last_frame(stacked_data))
end

function RLBase.reset!(p::StackFrames{T,N}) where {T,N}
    empty!(p.buffer)
    for _ in 1:capacity(p.buffer)
        push!(p.buffer, zeros(T, size(p.buffer)[1:N-1]))
    end
    p
end


text:
"""
    StaticInt(N::Int) -> StaticInt{N}()

A statically sized `Int`.
Use `StaticInt(N)` instead of `Val(N)` when you want it to behave like a number.
"""
struct StaticInt{N} <: Integer
    StaticInt{N}() where {N} = new{N::Int}()
end

const Zero = StaticInt{0}
const One = StaticInt{1}

Base.show(io::IO, ::StaticInt{N}) where {N} = print(io, "static($N)")

StaticInt(N::Int) = StaticInt{N}()
StaticInt(N::Integer) = StaticInt(convert(Int, N))
StaticInt(::StaticInt{N}) where {N} = StaticInt{N}()
StaticInt(::Val{N}) where {N} = StaticInt{N}()
# Base.Val(::StaticInt{N}) where {N} = Val{N}()
Base.convert(::Type{T}, ::StaticInt{N}) where {T<:Number,N} = convert(T, N)
Base.Bool(x::StaticInt{N}) where {N} = Bool(N)
Base.BigInt(x::StaticInt{N}) where {N} = BigInt(N)
Base.Integer(x::StaticInt{N}) where {N} = x
(::Type{T})(x::StaticInt{N}) where {T<:Integer,N} = T(N)
(::Type{T})(x::Int) where {T<:StaticInt} = StaticInt(x)
Base.convert(::Type{StaticInt{N}}, ::StaticInt{N}) where {N} = StaticInt{N}()

Base.promote_rule(::Type{<:StaticInt}, ::Type{T}) where {T<:Number} = promote_type(Int, T)
function Base.promote_rule(::Type{<:StaticInt}, ::Type{T}) where {T<:AbstractIrrational}
    return promote_type(Int, T)
end
# Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T <: AbstractIrrational} = promote_rule(T, Int)
for (S, T) in [(:Complex, :Real), (:Rational, :Integer), (:(Base.TwicePrecision), :Any)]
    @eval function Base.promote_rule(::Type{$S{T}}, ::Type{<:StaticInt}) where {T<:$T}
        return promote_type($S{T}, Int)
    end
end
function Base.promote_rule(::Type{Union{Nothing,Missing}}, ::Type{<:StaticInt})
    return Union{Nothing,Missing,Int}
end
function Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T>:Union{Missing,Nothing}}
    return promote_type(T, Int)
end
Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T>:Nothing} = promote_type(T, Int)
Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T>:Missing} = promote_type(T, Int)
for T in [:Bool, :Missing, :BigFloat, :BigInt, :Nothing, :Any]
    # let S = :Any
    @eval begin
        function Base.promote_rule(::Type{S}, ::Type{$T}) where {S<:StaticInt}
            return promote_type(Int, $T)
        end
        function Base.promote_rule(::Type{$T}, ::Type{S}) where {S<:StaticInt}
            return promote_type($T, Int)
        end
    end
end
Base.promote_rule(::Type{<:StaticInt}, ::Type{<:StaticInt}) = Int
Base.:(%)(::StaticInt{N}, ::Type{Integer}) where {N} = N

Base.eltype(::Type{T}) where {T<:StaticInt} = Int
Base.iszero(::Zero) = true
Base.iszero(::StaticInt) = false
Base.isone(::One) = true
Base.isone(::StaticInt) = false
Base.zero(::Type{T}) where {T<:StaticInt} = Zero()
Base.one(::Type{T}) where {T<:StaticInt} = One()

for T in [:Real, :Rational, :Integer]
    @eval begin
        @inline Base.:(+)(i::$T, ::Zero) = i
        @inline Base.:(+)(i::$T, ::StaticInt{M}) where {M} = i + M
        @inline Base.:(+)(::Zero, i::$T) = i
        @inline Base.:(+)(::StaticInt{M}, i::$T) where {M} = M + i
        @inline Base.:(-)(i::$T, ::Zero) = i
        @inline Base.:(-)(i::$T, ::StaticInt{M}) where {M} = i - M
        @inline Base.:(*)(i::$T, ::Zero) = Zero()
        @inline Base.:(*)(i::$T, ::One) = i
        @inline Base.:(*)(i::$T, ::StaticInt{M}) where {M} = i * M
        @inline Base.:(*)(::Zero, i::$T) = Zero()
        @inline Base.:(*)(::One, i::$T) = i
        @inline Base.:(*)(::StaticInt{M}, i::$T) where {M} = M * i
    end
end
@inline Base.:(+)(::Zero, ::Zero) = Zero()
@inline Base.:(+)(::Zero, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline Base.:(+)(::StaticInt{M}, ::Zero) where {M} = StaticInt{M}()

@inline Base.:(-)(::StaticInt{M}) where {M} = StaticInt{-M}()
@inline Base.:(-)(::StaticInt{M}, ::Zero) where {M} = StaticInt{M}()

@inline Base.:(*)(::Zero, ::Zero) = Zero()
@inline Base.:(*)(::One, ::Zero) = Zero()
@inline Base.:(*)(::Zero, ::One) = Zero()
@inline Base.:(*)(::One, ::One) = One()
@inline Base.:(*)(::StaticInt{M}, ::Zero) where {M} = Zero()
@inline Base.:(*)(::Zero, ::StaticInt{M}) where {M} = Zero()
@inline Base.:(*)(::StaticInt{M}, ::One) where {M} = StaticInt{M}()
@inline Base.:(*)(::One, ::StaticInt{M}) where {M} = StaticInt{M}()
for f in [:(+), :(-), :(*), :(/), :(÷), :(%), :(<<), :(>>), :(>>>), :(&), :(|), :(⊻)]
    @eval @generated function Base.$f(::StaticInt{M}, ::StaticInt{N}) where {M,N}
        return Expr(:call, Expr(:curly, :StaticInt, $f(M, N)))
    end
end
for f in [:(<<), :(>>), :(>>>)]
    @eval begin
        @inline Base.$f(::StaticInt{M}, x::UInt) where {M} = $f(M, x)
        @inline Base.$f(x::Integer, ::StaticInt{M}) where {M} = $f(x, M)
    end
end
for f in [:(==), :(!=), :(<), :(≤), :(>), :(≥)]
    @eval begin
        @inline Base.$f(::StaticInt{M}, ::StaticInt{N}) where {M,N} = $f(M, N)
        @inline Base.$f(::StaticInt{M}, x::Int) where {M} = $f(M, x)
        @inline Base.$f(x::Int, ::StaticInt{M}) where {M} = $f(x, M)
    end
end

@inline function maybe_static(f::F, g::G, x) where {F,G}
    L = f(x)
    if L === missing
        return g(x)
    else
        return static(L)
    end
end

@inline Base.widen(::StaticInt{N}) where {N} = widen(N)

Base.UnitRange{T}(start::StaticInt, stop) where {T<:Real} = UnitRange{T}(T(start), T(stop))
Base.UnitRange{T}(start, stop::StaticInt) where {T<:Real} = UnitRange{T}(T(start), T(stop))
function Base.UnitRange{T}(start::StaticInt, stop::StaticInt) where {T<:Real}
    return UnitRange{T}(T(start), T(stop))
end

Base.UnitRange(start::StaticInt, stop) = UnitRange(Int(start), stop)
Base.UnitRange(start, stop::StaticInt) = UnitRange(start, Int(stop))
Base.UnitRange(start::StaticInt, stop::StaticInt) = UnitRange(Int(start), Int(stop))

"""
    eq(x, y)

Equivalent to `!=` but if `x` and `y` are both static returns a `StaticBool.
"""
eq(x::X, y::Y) where {X,Y} = ifelse(is_static(X) & is_static(Y), static, identity)(x == y)
eq(x::X) where {X} = Base.Fix2(eq, x)

"""
    ne(x, y)

Equivalent to `!=` but if `x` and `y` are both static returns a `StaticBool.
"""
ne(x::X, y::Y) where {X,Y} = !eq(x, y)
ne(x::X) where {X} = Base.Fix2(ne, x)

"""
    gt(x, y)

Equivalent to `>` but if `x` and `y` are both static returns a `StaticBool.
"""
gt(x::X, y::Y) where {X,Y} = ifelse(is_static(X) & is_static(Y), static, identity)(x > y)
gt(x::X) where {X} = Base.Fix2(gt, x)

"""
    ge(x, y)

Equivalent to `>=` but if `x` and `y` are both static returns a `StaticBool.
"""
ge(x::X, y::Y) where {X,Y} = ifelse(is_static(X) & is_static(Y), static, identity)(x >= y)
ge(x::X) where {X} = Base.Fix2(ge, x)

"""
    le(x, y)

Equivalent to `<=` but if `x` and `y` are both static returns a `StaticBool.
"""
le(x::X, y::Y) where {X,Y} = ifelse(is_static(X) & is_static(Y), static, identity)(x <= y)
le(x::X) where {X} = Base.Fix2(le, x)

"""
    lt(x, y)

Equivalent to `<` but if `x` and `y` are both static returns a `StaticBool.
"""
lt(x::X, y::Y) where {X,Y} = ifelse(is_static(X) & is_static(Y), static, identity)(x < y)
lt(x::X) where {X} = Base.Fix2(lt, x)


text:
function fw_sim_01( semu_mass_ratios::Vector{Float64} ,  v::Vector{Float64} )
    #-----Init EMRs-------------------------------------------------------------------------------------
    mr_sA_2_M_0 = semu_mass_ratios[1]
    mr_sA_2_M_1 = semu_mass_ratios[2]
    mr_sA_3_M_0 = semu_mass_ratios[3]
    mr_sA_3_M_1 = semu_mass_ratios[4]
    mr_sA_2_3_M_0 = semu_mass_ratios[5]
    mr_sA_2_3_M_1 = semu_mass_ratios[6]
    mr_sA_2_3_M_2 = semu_mass_ratios[7]
    mr_sA_1_2_3_M_0 = semu_mass_ratios[8]
    mr_sA_1_2_3_M_1 = semu_mass_ratios[9]
    mr_sA_1_2_3_M_2 = semu_mass_ratios[10]
    mr_sA_1_2_3_M_3 = semu_mass_ratios[11]
    #-----Init Fluxes-----------------------------------------------------------------------------------
    v_R01 = v[1]
    v_R02 = v[2]
    v_R03 = v[3]
    v_R04 = v[4]
    v_R05 = v[5]
    v_R06 = v[6]
    #---------------------------------------------------------------------------------------------------
    #----Init A---------------------------------------------------------------------------------
    sA_1_1::Array{Float64,2} = zeros(Float64,5,5)
    sA_1_1[1,1] = ( (  -1.0  * v_R01   )  +(  -1.0  * v_R03   )   )
    sA_1_1[3,1] =  v_R04
    sA_1_1[4,1] =  v_R02
    sA_1_1[2,2] = ( (  -1.0  * v_R01   )  +(  -1.0  * v_R03   )   )
    sA_1_1[4,2] =  v_R05
    sA_1_1[5,2] =  v_R02
    sA_1_1[3,3] = ( (  -1.0  * v_R04   )   )
    sA_1_1[5,3] =  v_R05
    sA_1_1[1,4] =  v_R03
    sA_1_1[4,4] = ( (  -1.0  * v_R02   )  +(  -1.0  * v_R05   )   )
    sA_1_1[2,5] =  v_R03
    sA_1_1[5,5] = ( (  -1.0  * v_R05   )  +(  -1.0  * v_R02   )   )
    #----Init B---------------------------------------------------------------------------
    sB_1_1::Array{Float64,2} =  zeros(Float64,5,2)
    sB_1_1[2,1] = (  -1.0  * v_R01   )
    sB_1_1[1,2] = (  -1.0  * v_R01   )
    #----Init Y---------------------------------------------------------------------------
    sY_1_1::Array{Float64,2} =  zeros(Float64,2,2)
    sY_1_1[1,1] = ( (  mr_sA_3_M_0   )   )
    sY_1_1[2,1] = ( (  mr_sA_2_M_0   )   )
    sY_1_1[1,2] = ( (  mr_sA_3_M_1   )   )
    sY_1_1[2,2] = ( (  mr_sA_2_M_1   )   )

    #------COMPUTE!-----------------------------------------------------------------------
    sX_1_1 = sA_1_1 \ ( sB_1_1 * sY_1_1 )
    #-------------------------------------------------------------------------------------
    mr_sB_2_M_0 = sX_1_1[1,1]
    mr_sB_3_M_0 = sX_1_1[2,1]
    mr_sC_1_M_0 = sX_1_1[3,1]
    mr_sD_2_M_0 = sX_1_1[4,1]
    mr_sD_3_M_0 = sX_1_1[5,1]
    mr_sB_2_M_1 = sX_1_1[1,2]
    mr_sB_3_M_1 = sX_1_1[2,2]
    mr_sC_1_M_1 = sX_1_1[3,2]
    mr_sD_2_M_1 = sX_1_1[4,2]
    mr_sD_3_M_1 = sX_1_1[5,2]
    #-------------------------------------------------------------------------------------
    #----Init A---------------------------------------------------------------------------------
    sA_2_1::Array{Float64,2} =  zeros(Float64,2,2)
    sA_2_1[1,1] = ( (  -1.0  * v_R01   )  +(  -1.0  * v_R03   )   )
    sA_2_1[2,1] =  v_R02
    sA_2_1[1,2] =  v_R03
    sA_2_1[2,2] = ( (  -1.0  * v_R02   )  +(  -1.0  * v_R05   )   )
    #----Init B---------------------------------------------------------------------------
    sB_2_1::Array{Float64,2} =  zeros(Float64,2,2)
    sB_2_1[1,1] = (  -1.0  * v_R01   )
    sB_2_1[2,2] = (  -1.0  * v_R05   )
    #----Init Y---------------------------------------------------------------------------
    sY_2_1::Array{Float64,2} =  zeros(Float64,2,3)
    sY_2_1[1,1] = ( (  mr_sA_2_3_M_0   )   )
    sY_2_1[2,1] = ( (  mr_sB_3_M_0  * mr_sC_1_M_0   )   )
    sY_2_1[1,2] = ( (  mr_sA_2_3_M_1   )   )
    sY_2_1[2,2] = ( (  mr_sB_3_M_0  * mr_sC_1_M_1   )  +(  mr_sB_3_M_1  * mr_sC_1_M_0   )   )
    sY_2_1[1,3] = ( (  mr_sA_2_3_M_2   )   )
    sY_2_1[2,3] = ( (  mr_sB_3_M_1  * mr_sC_1_M_1   )   )

    #------COMPUTE!-----------------------------------------------------------------------
    sX_2_1 = sA_2_1 \ ( sB_2_1 * sY_2_1 )
    #-------------------------------------------------------------------------------------
    mr_sB_2_3_M_0 = sX_2_1[1,1]
    mr_sD_2_3_M_0 = sX_2_1[2,1]
    mr_sB_2_3_M_1 = sX_2_1[1,2]
    mr_sD_2_3_M_1 = sX_2_1[2,2]
    mr_sB_2_3_M_2 = sX_2_1[1,3]
    mr_sD_2_3_M_2 = sX_2_1[2,3]
    #-------------------------------------------------------------------------------------
    #----Init A---------------------------------------------------------------------------------
    sA_3_1::Array{Float64,2} =  zeros(Float64,3,3)
    sA_3_1[1,1] = ( (  -1.0  * v_R03   )  +(  -1.0  * v_R01   )   )
    sA_3_1[2,1] =  v_R02
    sA_3_1[1,2] =  v_R03
    sA_3_1[2,2] = ( (  -1.0  * v_R05   )  +(  -1.0  * v_R02   )   )
    sA_3_1[3,2] =  v_R06
    sA_3_1[3,3] = ( (  -1.0  * v_R06   )   )
    #----Init B---------------------------------------------------------------------------
    sB_3_1::Array{Float64,2} =  zeros(Float64,3,2)
    sB_3_1[2,1] = (  -1.0  * v_R05   )
    sB_3_1[1,2] = (  -1.0  * v_R01   )
    #----Init Y---------------------------------------------------------------------------
    sY_3_1::Array{Float64,2} =  zeros(Float64,2,4)
    sY_3_1[1,1] = ( (  mr_sB_2_3_M_0  * mr_sC_1_M_0   )   )
    sY_3_1[2,1] = ( (  mr_sA_1_2_3_M_0   )   )
    sY_3_1[1,2] = ( (  mr_sB_2_3_M_0  * mr_sC_1_M_1   )  +(  mr_sB_2_3_M_1  * mr_sC_1_M_0   )   )
    sY_3_1[2,2] = ( (  mr_sA_1_2_3_M_1   )   )
    sY_3_1[1,3] = ( (  mr_sB_2_3_M_1  * mr_sC_1_M_1   )  +(  mr_sB_2_3_M_2  * mr_sC_1_M_0   )   )
    sY_3_1[2,3] = ( (  mr_sA_1_2_3_M_2   )   )
    sY_3_1[1,4] = ( (  mr_sB_2_3_M_2  * mr_sC_1_M_1   )   )
    sY_3_1[2,4] = ( (  mr_sA_1_2_3_M_3   )   )

    #------COMPUTE!-----------------------------------------------------------------------
    sX_3_1 = sA_3_1 \ ( sB_3_1 * sY_3_1 )
    #-------------------------------------------------------------------------------------
    mr_sB_1_2_3_M_0 = sX_3_1[1,1]
    mr_sD_1_2_3_M_0 = sX_3_1[2,1]
    mr_sF_1_2_3_M_0 = sX_3_1[3,1]
    mr_sB_1_2_3_M_1 = sX_3_1[1,2]
    mr_sD_1_2_3_M_1 = sX_3_1[2,2]
    mr_sF_1_2_3_M_1 = sX_3_1[3,2]
    mr_sB_1_2_3_M_2 = sX_3_1[1,3]
    mr_sD_1_2_3_M_2 = sX_3_1[2,3]
    mr_sF_1_2_3_M_2 = sX_3_1[3,3]
    mr_sB_1_2_3_M_3 = sX_3_1[1,4]
    mr_sD_1_2_3_M_3 = sX_3_1[2,4]
    mr_sF_1_2_3_M_3 = sX_3_1[3,4]
    #-------------------------------------------------------------------------------------
    #---RESULTS-----------------------------------------------------------------------------------------
    target_smrs = Array{Float64,1}(4)
    target_smrs[1] = mr_sF_1_2_3_M_0
    target_smrs[2] = mr_sF_1_2_3_M_1
    target_smrs[3] = mr_sF_1_2_3_M_2
    target_smrs[4] = mr_sF_1_2_3_M_3
    #---------------------------------------------------------------------------------------------------
    return target_smrs
    #---------------------------------------------------------------------------------------------------
end

semu_mass_ratios = [0.;1;1;0;0;1;0;0;1;0;0]

#                    v1  v2     v3   v4   v5   v6      Ain  Eout  Fout
reference_flux = [  100.; 110 ; 50 ; 20 ; 20 ; 80   ;  100 ; 60 ; 80 ]

fw_result = fw_sim_01(semu_mass_ratios,reference_flux)


text:
# Tests for data project functionality

using Pkg

using DataSets:
    TomlFileDataProject,
    ActiveDataProject,
    StackedDataProject,
    project_name

test_project_names = ["a_text_file",
                      "a_tree_example",
                      "embedded_blob",
                      "embedded_tree",
                      "gzipped_table",
                      "old_backend_blob",
                      "old_backend_tree",
                      "some_namespace/a_text_file"]

@testset "TomlFileDataProject" begin
    proj = TomlFileDataProject(abspath("Data.toml"))
    # getindex, get
    @test proj["a_text_file"].uuid == UUID("b498f769-a7f6-4f67-8d74-40b770398f26")
    @test_throws ErrorException proj["nonexistent_data"]
    @test isnothing(get(proj, "nonexistent_data", nothing))

    # keys
    @test sort(collect(keys(proj))) == test_project_names
    @test haskey(proj, "a_text_file")
    @test !haskey(proj, "nonexistent_data")

    # iteration
    @test sort(getproperty.(collect(proj), :name)) == test_project_names
    @test sort(first.(pairs(proj))) == test_project_names

    # identity
    @test project_name(proj) == abspath("Data.toml")

    # Test @__DIR__ templating
    # Use `cleanpath` as there's currently a mixture of / and \ on windows
    # which does work, but is quite ugly.
    cleanpath(p) = replace(p, '\\'=>'/')
    @test cleanpath(proj["a_text_file"].storage["path"]) == cleanpath(joinpath(@__DIR__, "data", "file.txt"))
end

@testset "TomlFileDataProject live updates" begin
    # Test live updating when the file is rewritten
    mktemp() do path,io
        write(io, """
        data_config_version=0

        [[datasets]]
        description="A text file"
        name="a_text_file"
        uuid="b498f769-a7f6-4f67-8d74-40b770398f26"

            [datasets.storage]
            driver="FileSystem"
            type="Blob"
            path="@__DIR__/data/file.txt"
        """)
        flush(io)

        proj = TomlFileDataProject(path)

        @test collect(keys(proj)) == ["a_text_file"]
        @test proj["a_text_file"].uuid == UUID("b498f769-a7f6-4f67-8d74-40b770398f26")

        write(io, """

        [[datasets]]
        description="A text file 2"
        name="a_text_file_2"
        uuid="58992dd5-86ed-4747-b9fc-320e13a03504"

            [datasets.storage]
            driver="FileSystem"
            type="Blob"
            path="@__DIR__/data/file2.txt"
        """)
        flush(io)

        @test sort(collect(keys(proj))) == ["a_text_file", "a_text_file_2"]
        @test proj["a_text_file_2"].uuid == UUID("58992dd5-86ed-4747-b9fc-320e13a03504")
    end
end

@testset "ActiveDataProject" begin
    proj = ActiveDataProject()
    @test project_name(proj) == nothing
    @test collect(keys(proj)) == []
    @test collect(pairs(proj)) == []
    @test collect(proj) == []
    proj_dir = joinpath(@__DIR__, "active_project")
    Pkg.activate(proj_dir)
    @test project_name(proj) == joinpath(proj_dir, "Data.toml")
    @test collect(keys(proj)) == ["a_text_file"]
    @test proj["a_text_file"].uuid == UUID("314996ef-12be-40d0-912c-9755af354fdb")
    Pkg.activate()
    @test project_name(proj) == nothing
    @test collect(keys(proj)) == []
end

@testset "StackedDataProject" begin
    proj = StackedDataProject()
    @test collect(keys(proj)) == []

    push!(proj, TomlFileDataProject(joinpath(@__DIR__, "active_project", "Data.toml")))
    push!(proj, TomlFileDataProject(joinpath(@__DIR__, "Data.toml")))

    @test sort(collect(keys(proj))) == test_project_names
    # Data "a_text_file" should be found in the first project in the stack,
    # overriding the data of the same name in the second project.
    @test proj["a_text_file"].uuid == UUID("314996ef-12be-40d0-912c-9755af354fdb")
end

@testset "JULIA_DATASETS_PATH" begin
    # Test JULIA_DATASETS_PATH environment variable and path unpacking
    paths_sep = Sys.iswindows() ? ';' : ':'
    datasets_paths = join(["@",
                           joinpath(@__DIR__, "Data.toml"),
                           ""], paths_sep)
    fake_env = Dict("JULIA_DATASETS_PATH"=>datasets_paths)
    proj = DataSets.create_project_stack(fake_env)
    @test proj.projects[1] isa ActiveDataProject

    @test proj.projects[2] isa TomlFileDataProject
    @test project_name(proj.projects[2]) == joinpath(@__DIR__, "Data.toml")

    @test proj.projects[3] isa TomlFileDataProject
    @test project_name(proj.projects[3]) == joinpath(homedir(), ".julia", "datasets", "Data.toml")

    # Test that __init__ takes global DataSets.PROJECT from ENV
    empty!(DataSets.PROJECT)
    ENV["JULIA_DATASETS_PATH"] = @__DIR__
    DataSets.__init__()
    @test DataSets.PROJECT.projects[1] isa TomlFileDataProject
    @test project_name(DataSets.PROJECT.projects[1]) == joinpath(@__DIR__, "Data.toml")
end


text:
#+Tag L:target
# This is an variable example [[target][]]
A=1


text:
"""
```julia
three = 3
3
markable_three = Unmarked(three)
3
isunmarked(markable_three)
true
@mark!(markable_three)
3
ismarked(markable_three)
true
```
"""
macro mark!(x)
    quote
        $(esc(x)) = reinterpret(typeof($(esc(x))), lsbit(typeof($(esc(x)))) | ($(esc(x))))
    end
end

"""
```julia
three = 3
3
markable_three = Unmarked(three)
3
isunmarked(markable_three)
true
@mark!(markable_three)
3
ismarked(markable_three)
true
@unmark!(markable_three)
3
ismarked(markable_three)
false
```
"""
macro unmark!(x)
    quote
        $(esc(x)) = reinterpret(typeof($(esc(x))), msbitsof($(esc(x))))
    end
end

for (M,I) in ((:MarkInt128, :Int128), (:MarkInt64, :Int64),
              (:MarkInt32, :Int32), (:MarkInt16, :Int16), (:MarkInt8, :Int8),
              (:MarkUInt128, :UInt128), (:MarkUInt64, :UInt64),
              (:MarkUInt32, :UInt32), (:MarkUInt16, :UInt16), (:MarkUInt8, :UInt8))
  @eval begin
    @inline ismarked(x::$M) = isodd(reinterpret($I, x))
    @inline isunmarked(x::$M) = iseven(reinterpret($I, x))
    @inline ismarked(x::$I) = false
    @inline isunmarked(x::$I) = true

    @inline mtype(::Type{$I}) = $M
    @inline mtype(::Type{$M}) = $M
    @inline itype(::Type{$I}) = $I
    @inline itype(::Type{$M}) = $I
    @inline mtype(x::$I) = $M
    @inline mtype(x::$M) = $M
    @inline itype(x::$I) = $I
    @inline itype(x::$M) = $I
    @inline mtyped(x::$I) = reinterpret($M, x)
    @inline mtyped(x::$M) = x
    @inline ityped(x::$M) = reinterpret($I, x)
    @inline ityped(x::$I) = x

    @inline mark_unmarked(x::$M) = reinterpret($M, (reinterpret($I, x) | lsbit($I)))
    @inline unmark_marked(x::$M) = reinterpret($M, (msbitsof(reinterpret($I, x))))
    @inline mark(x::$M) = ismarked(x) ? x : mark_unmarked(x)
    @inline unmark(x::$M) = isunmarked(x) ? x : unmark_marked(x)
  end
end

find_marked(x::Vector{T}) where {T<:MarkableInteger}  = findall(map(ismarked, x))
find_unmarked(x::Vector{T}) where {T<:MarkableInteger}  = findall(map(isunmarked, x))
all_marked(x::Vector{T}) where {T<:MarkableInteger} = x[find_marked(x)]
all_unmarked(x::Vector{T}) where {T<:MarkableInteger} = x[find_unmarked(x)]