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

################################### Unit Commitment tests ##################################
@testset "Thermal UC With DC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
    uc_constraint_names = [:ramp_up_ThermalStandard,
                           :ramp_dn_ThermalStandard,
                           :duration_up_ThermalStandard,
                           :duration_dn_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment)

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

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

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

@testset "Thermal UC With AC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
    uc_constraint_names = [:ramp_up_ThermalStandard,
                           :ramp_dn_ThermalStandard,
                           :duration_up_ThermalStandard,
                           :duration_dn_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment)

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

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

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

################################### Basic Unit Commitment tests ############################
@testset "Thermal Basic UC With DC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalBasicUnitCommitment)

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

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

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

@testset "Thermal Basic UC With AC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
     model = DeviceModel(PSY.ThermalStandard, PSI.ThermalBasicUnitCommitment)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

############################# UC validation tests ##########################################
branches = Dict{Symbol, PSI.DeviceModel}()
services = Dict{Symbol, PSI.ServiceModel}()
ED_devices = Dict{Symbol, DeviceModel}(:Generators => PSI.DeviceModel(PSY.ThermalStandard, PSI.ThermalRampLimited),
                                        :Loads =>  PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
UC_devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment),
                                        :Loads =>  DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
# Testing Ramping Constraint
@testset "Solving UC with CopperPlate for testing Ramping Constraints" begin
node = Bus(1,"nodeA", "PV", 0, 1.0, (min = 0.9, max=1.05), 230)
load = PowerLoad("Bus1", true, node,nothing, 0.4, 0.9861, 1.0, 2.0)
    DA_ramp = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  4:00:00", "d/m/y  H:M:S"))
    gen_ramp = [ThermalStandard("Alta", true, node,0.20, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min=0.0, max=0.40), nothing,
                            nothing, nothing),
            ThreePartCost((0.0, 1400.0), 0.0, 4.0, 2.0)
            ),
            ThermalStandard("Park City", true, node,0.70,0.20,
                TechThermal(2.0,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min=0.7, max=2.20), nothing,
                            (up = 0.010625, down= 0.010625),nothing),
                ThreePartCost((0.0, 1500.0), 0.0, 1.5, 0.75)
            )];
    ramp_load = [ 0.9, 1.1, 2.485, 2.175, 0.9];
    load_forecast_ramp = Deterministic("get_maxactivepower", TimeArray(DA_ramp, ramp_load))
    ramp_test_sys = PSY.System(100.0)
    add_component!(ramp_test_sys, node)
    add_component!(ramp_test_sys, load)
    add_component!(ramp_test_sys, gen_ramp[1])
    add_component!(ramp_test_sys, gen_ramp[2])
    add_forecast!(ramp_test_sys, load, load_forecast_ramp)

    template = OperationsTemplate(CopperPlatePowerModel, ED_devices, branches, services)
    ED = OperationsProblem(TestOpProblem, template,
                        ramp_test_sys; optimizer = Cbc_optimizer,
                        use_parameters = true)
    psi_checksolve_test(ED, [MOI.OPTIMAL], 11191.00)
    moi_tests(ED, true, 10, 10, 10, 0, 5, false)
end

# Testing Duration Constraints
@testset "Solving UC with CopperPlate for testing Duration Constraints" begin
node = Bus(1,"nodeA", "PV", 0, 1.0, (min = 0.9, max=1.05), 230)
load = PowerLoad("Bus1", true, node,nothing, 0.4, 0.9861, 1.0, 2.0)
    DA_dur  = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  6:00:00", "d/m/y  H:M:S"))
    gens_dur = [ThermalStandard("Alta", true, node,0.40, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min=0.3, max=0.9), nothing,
                            nothing, (up=4, down=2)),
            ThreePartCost((0.0, 1400.0), 0.0, 4.0, 2.0)
            ),
            ThermalStandard("Park City", true, node,1.70,0.20,
                TechThermal(2.2125,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                                (min=0.7, max=2.2), nothing,
                                nothing, (up=6, down=4)),
                ThreePartCost((0.0, 1500.0), 0.0, 1.5, 0.75)
            )];

    duration_load = [0.3, 0.6, 0.8, 0.7, 1.7, 0.9, 0.7]
    load_forecast_dur = Deterministic("get_maxactivepower", TimeArray(DA_dur, duration_load))
    duration_test_sys = PSY.System(100.0)
    add_component!(duration_test_sys, node)
    add_component!(duration_test_sys, load)
    add_component!(duration_test_sys, gens_dur[1])
    add_component!(duration_test_sys, gens_dur[2])
    add_forecast!(duration_test_sys, load, load_forecast_dur)

    status = [1.0,0.0]
    up_time = [2.0,0.0]
    down_time = [0.0,3.0]

    alta = gens_dur[1]
    init_cond = PSI.DICKDA()
    init_cond[PSI.ICKey(PSI.DeviceStatus,typeof(alta))] = build_init(gens_dur, status)
    init_cond[PSI.ICKey(PSI.TimeDurationON,typeof(alta))] = build_init(gens_dur, up_time)
    init_cond[PSI.ICKey(PSI.TimeDurationOFF,typeof(alta))] = build_init(gens_dur, down_time)

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

## PWL linear Cost implementation test
@testset "Solving UC with CopperPlate testing Linear PWL" begin
node = Bus(1,"nodeA", "PV", 0, 1.0, (min = 0.9, max=1.05), 230)
load = PowerLoad("Bus1", true, node,nothing, 0.4, 0.9861, 1.0, 2.0)
    gens_cost = [ThermalStandard("Alta", true, node,0.52, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.22, max = 0.55), nothing,
                            nothing, nothing),
            ThreePartCost([ (589.99, 0.220),(884.99, 0.33)
                    ,(1210.04, 0.44),(1543.44, 0.55)],532.44, 5665.23, 0.0)
            ),
            ThermalStandard("Park City", true, node,0.62,0.20,
                TechThermal(2.2125,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.62, max = 1.55), nothing,
                            nothing, nothing),
                ThreePartCost([   (1264.80, 0.62),(1897.20, 0.93),
                    (2594.4787, 1.24),(3433.04, 1.55)   ], 235.397, 5665.23, 0.0)
            )];
    DA_cost  = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  1:00:00", "d/m/y  H:M:S"))
    cost_load = [1.3,2.1];
    load_forecast_cost = Deterministic("get_maxactivepower", TimeArray(DA_cost, cost_load))
    cost_test_sys = PSY.System(100.0)
    add_component!(cost_test_sys, node)
    add_component!(cost_test_sys, load)
    add_component!(cost_test_sys, gens_cost[1])
    add_component!(cost_test_sys, gens_cost[2])
    add_forecast!(cost_test_sys, load, load_forecast_cost)

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

## PWL SOS-2 Cost implementation test
@testset "Solving UC with CopperPlate testing SOS2 implementation" begin
    node = Bus(1,"nodeA", "PV", 0, 1.0, (min = 0.9, max=1.05), 230)
    load = PowerLoad("Bus1", true, node,nothing, 0.4, 0.9861, 1.0, 2.0)
    gens_cost_sos = [ThermalStandard("Alta", true, node,0.52, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.22, max = 0.55), nothing,
                            nothing, nothing),
            ThreePartCost([ (1122.43, 0.22),(1417.43, 0.33),
                    (1742.48, 0.44),(2075.88, 0.55) ],0.0, 5665.23, 0.0)
            ),
            ThermalStandard("Park City", true, node,0.62,0.20,
                TechThermal(2.2125,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.62, max = 1.55), nothing,
                            nothing, nothing),
                ThreePartCost([ (1500.19, 0.62),(2132.59, 0.929),
                    (2829.875, 1.24),(3668.444, 1.55)], 0.0, 5665.23, 0.0)
            )];
    DA_cost_sos   = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  1:00:00", "d/m/y  H:M:S"))
    cost_sos_load = [1.3,2.1];
    load_forecast_cost_sos  = Deterministic("get_maxactivepower", TimeArray(DA_cost_sos, cost_sos_load))
    cost_test_sos_sys = PSY.System(100.0)
    add_component!(cost_test_sos_sys, node)
    add_component!(cost_test_sos_sys, load)
    add_component!(cost_test_sos_sys, gens_cost_sos[1])
    add_component!(cost_test_sos_sys, gens_cost_sos[2])
    add_forecast!(cost_test_sos_sys, load, load_forecast_cost_sos)

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


text:
using POMDPGym, POMDPs, Test
mdp = CartPoleMDP()

s = rand(initialstate(mdp))
@test length(s) == 5
@test s[end] == 0f0

sp, r = gen(mdp, s, -1)

ext: jl

max_forks_repo_name: ancorso/POMDPGym


text:
using Cthulhu, Test, DeepDiffs

include("utils.jl")

m = Module()
Core.eval(m, Meta.parseall(raw"""
function foo(x, y)
    z = x + y
    if z < 4
        z += 1
    end
    u = (x -> x + z)(x)
    v = Ref{Union{Int, Missing}}(x)[] + y
    return u + v
end
"""; filename="foobar.jl"))

@testset "IRShow tests" begin
    tf = (true, false)

    @testset "optimize: $optimize" for optimize in tf
        _, src, infos, mi, rt, slottypes = process(m.foo, (Int, Int); optimize);

        @testset "debuginfo: $debuginfo" for debuginfo in instances(Cthulhu.DInfo.DebugInfo)
            @testset "iswarn: $iswarn" for iswarn in tf
                @testset "hide_type_stable: $hide_type_stable" for hide_type_stable in tf
                    @testset "inline_cost: $inline_cost" for inline_cost in tf
                        @testset "type_annotations: $type_annotations" for type_annotations in tf
                            !optimize && debuginfo === Cthulhu.DInfo.compact && continue
                            !optimize && inline_cost && continue

                            s = sprint(; context=:color=>true) do io
                                Cthulhu.cthulhu_typed(io, debuginfo,
                                                      src, rt, mi;
                                                      iswarn, hide_type_stable, inline_cost, type_annotations)
                            end
                            s = strip_base_linenums(s)

                            ground_truth = read(irshow_filename(optimize, debuginfo, iswarn, hide_type_stable, inline_cost, type_annotations), String)
                            if Sys.iswindows()
                                ground_truth = replace(ground_truth, "\r\n" => "\n")
                            end
                            @test s == ground_truth
                            s != ground_truth && println(deepdiff(s, ground_truth))
                        end
                    end
                end
            end
        end
    end
end


text:
# Decl
function dumpColor(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_dumpColor(x)
end

function getBeginLoc(x::AbstractDecl)
    @check_ptrs x
    return SourceLocation(clang_Decl_getBeginLoc(x))
end

function getEndLoc(x::AbstractDecl)
    @check_ptrs x
    return SourceLocation(clang_Decl_getEndLoc(x))
end

function getASTContext(x::AbstractDecl)
    @check_ptrs x
    return ASTContext(clang_Decl_getASTContext(x))
end

function isOutOfLine(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isOutOfLine(x)
end

function getLocation(x::AbstractDecl)
    @check_ptrs x
    return SourceLocation(clang_Decl_getLocation(x))
end

function getDeclKindName(x::AbstractDecl)
    @check_ptrs x
    return unsafe_string(clang_Decl_getDeclKindName(x))
end

function getNextDeclInContext(x::AbstractDecl)
    @check_ptrs x
    return Decl(clang_Decl_getNextDeclInContext(x))
end

function getDeclContext(x::AbstractDecl)
    @check_ptrs x
    return DeclContext(clang_Decl_getDeclContext(x))
end

function getNonClosureContext(x::AbstractDecl)
    @check_ptrs x
    return Decl(clang_Decl_getNonClosureContext(x))
end

function isInAnonymousNamespace(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isInAnonymousNamespace(x)
end

function isInStdNamespace(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isInStdNamespace(x)
end

function getLangOpts(x::AbstractDecl)
    @check_ptrs x
    return LangOptions(clang_Decl_getLangOpts(x))
end

function getLexicalDeclContext(x::AbstractDecl)
    @check_ptrs x
    return DeclContext(clang_Decl_getLexicalDeclContext(x))
end

function isTemplated(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isTemplated(x)
end

function isCanonicalDecl(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isCanonicalDecl(x)
end

function getPreviousDecl(x::AbstractDecl)
    @check_ptrs x
    return Decl(clang_Decl_getPreviousDecl(x))
end

function isFirstDecl(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isFirstDecl(x)
end

function getMostRecentDecl(x::AbstractDecl)
    @check_ptrs x
    return Decl(clang_Decl_getMostRecentDecl(x))
end

function isTemplateParameter(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isTemplateParameter(x)
end

function isTemplateDecl(x::AbstractDecl)
    @check_ptrs x
    return clang_Decl_isTemplateDecl(x)
end

function getDescribedTemplate(x::AbstractDecl)
    @check_ptrs x
    return TemplateDecl(clang_Decl_getDescribedTemplate(x))
end

function setDeclContext(decl::AbstractDecl, x::DeclContext)
    @check_ptrs x decl
    return clang_Decl_setDeclContext(decl, x)
end

function setLexicalDeclContext(decl::AbstractDecl, x::DeclContext)
    @check_ptrs x decl
    return clang_Decl_setLexicalDeclContext(decl, x)
end

function ClassTemplateDecl(x::AbstractDecl)
    @check_ptrs x
    return ClassTemplateDecl(clang_Decl_castToClassTemplateDecl(x))
end

function ValueDecl(x::AbstractDecl)
    @check_ptrs x
    return ValueDecl(clang_Decl_castToValueDecl(x))
end

# NamedDecl
function isOutOfLine(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_isOutOfLine(x)
end

function getIdentifier(x::AbstractNamedDecl)
    @check_ptrs x
    return IdentifierInfo(clang_NamedDecl_getIdentifier(x))
end

function getName(x::AbstractNamedDecl)
    @check_ptrs x
    return unsafe_string(clang_NamedDecl_getName(x))
end

function getDeclName(x::AbstractNamedDecl)
    @check_ptrs x
    return DeclarationName(clang_NamedDecl_getDeclName(x))
end

function setDeclName(x::AbstractNamedDecl, name::DeclarationName)
    @check_ptrs x
    return clang_NamedDecl_setDeclName(x, name)
end

function hasLinkage(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_hasLinkage(x)
end

function isCXXClassMember(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_isCXXClassMember(x)
end

function isCXXInstanceMember(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_isCXXInstanceMember(x)
end

function hasExternalFormalLinkage(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_hasExternalFormalLinkage(x)
end

function isExternallyVisible(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_isExternallyVisible(x)
end

function isExternallyDeclarable(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_isExternallyDeclarable(x)
end

function getUnderlyingDecl(x::AbstractNamedDecl)
    @check_ptrs x
    return NamedDecl(clang_NamedDecl_getUnderlyingDecl(x))
end

function getMostRecentDecl(x::AbstractNamedDecl)
    @check_ptrs x
    return clang_NamedDecl_getMostRecentDecl(x)
end

function TypeDecl(x::NamedDecl)
    @check_ptrs x
    return TypeDecl(clang_NamedDecl_castToTypeDecl(x))
end

# ValueDecl
function getType(x::AbstractValueDecl)
    @check_ptrs x
    return QualType(clang_ValueDecl_getType(x))
end

function setType(x::AbstractValueDecl, ty::AbstractQualType)
    @check_ptrs x
    return clang_ValueDecl_setType(x, ty)
end

function isWeak(x::AbstractValueDecl)
    @check_ptrs x
    return clang_ValueDecl_isWeak(x)
end

# TypeDecl
function getTypeForDecl(x::AbstractTypeDecl)::CXType_
    @check_ptrs x
    return clang_TypeDecl_getTypeForDecl(x)
end
getTypeForDecl(x::NamedDecl) = getTypeForDecl(TypeDecl(x))

function setTypeForDecl(x::AbstractTypeDecl, ty_ptr::CXType_)
    @check_ptrs x
    return clang_TypeDecl_setTypeForDecl(x, ty_ptr)
end
setTypeForDecl(x::AbstractTypeDecl, ty::AbstractQualType) = setTypeForDecl(x, get_type_ptr(ty))

function getBeginLoc(x::AbstractTypeDecl)
    @check_ptrs x
    return SourceLocation(clang_TypeDecl_getBeginLoc(x))
end

function setLocStart(x::AbstractTypeDecl, loc::SourceLocation)
    @check_ptrs x
    return clang_TypeDecl_setLocStart(x, loc)
end

# TypedefNameDecl
function getUnderlyingType(x::AbstractTypedefNameDecl)
    @check_ptrs x
    return QualType(clang_TypedefNameDecl_getUnderlyingType(x))
end

function getCanonicalDecl(x::AbstractTypedefNameDecl)
    @check_ptrs x
    return TypedefNameDecl(clang_TypedefNameDecl_getCanonicalDecl(x))
end

function getAnonDeclWithTypedefName(x::AbstractTypedefNameDecl, any_redecl::Bool=false)
    @check_ptrs x
    return TagDecl(clang_TypedefNameDecl_getAnonDeclWithTypedefName(x, any_redecl))
end

function isTransparentTag(x::AbstractTypedefNameDecl)
    @check_ptrs x
    return clang_TypedefNameDecl_isTransparentTag(x)
end

# TagDecl
function DeclContext(x::AbstractTagDecl)
    @check_ptrs x
    return DeclContext(clang_TagDecl_castToDeclContext(x))
end

function getCanonicalDecl(x::AbstractTagDecl)
    @check_ptrs x
    return TagDecl(clang_TagDecl_getCanonicalDecl(x))
end

function isThisDeclarationADefinition(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isThisDeclarationADefinition(x)
end

function isCompleteDefinition(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isCompleteDefinition(x)
end

function setCompleteDefinition(x::AbstractTagDecl, v::Bool)
    @check_ptrs x
    return clang_TagDecl_setCompleteDefinition(x, v)
end

function isBeingDefined(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isBeingDefined(x)
end

function isFreeStanding(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isFreeStanding(x)
end

function startDefinition(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_startDefinition(x)
end

function getDefinition(x::AbstractTagDecl)
    @check_ptrs x
    return TagDecl(clang_TagDecl_getDefinition(x))
end

function getKindName(x::AbstractTagDecl)
    @check_ptrs x
    return unsafe_string(clang_TagDecl_getKindName(x))
end

function getTagKind(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_getTagKind(x)
end

function isStruct(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isStruct(x)
end

function isInterface(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isInterface(x)
end

function isClass(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isClass(x)
end

function isUnion(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isUnion(x)
end

function isEnum(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_isEnum(x)
end

function hasNameForLinkage(x::AbstractTagDecl)
    @check_ptrs x
    return clang_TagDecl_hasNameForLinkage(x)
end

function getTypedefNameForAnonDecl(x::AbstractTagDecl)
    @check_ptrs x
    return TypedefNameDecl(clang_TagDecl_getTypedefNameForAnonDecl(x))
end

function getQualifier(x::AbstractTagDecl)
    @check_ptrs x
    return NestedNameSpecifier(clang_TagDecl_getQualifier(x))
end

# EnumDecl
function getCanonicalDecl(x::EnumDecl)
    @check_ptrs x
    return EnumDecl(clang_EnumDecl_getCanonicalDecl(x))
end

function getPreviousDecl(x::EnumDecl)
    @check_ptrs x
    return EnumDecl(clang_EnumDecl_getPreviousDecl(x))
end

function getMostRecentDecl(x::EnumDecl)
    @check_ptrs x
    return EnumDecl(clang_EnumDecl_getMostRecentDecl(x))
end

function getDefinition(x::EnumDecl)
    @check_ptrs x
    return EnumDecl(clang_EnumDecl_getDefinition(x))
end

function getIntegerType(x::EnumDecl)
    @check_ptrs x
    return QualType(clang_EnumDecl_getIntegerType(x))
end

# RecordDecl
function getPreviousDecl(x::AbstractRecordDecl)
    @check_ptrs x
    return RecordDecl(clang_RecordDecl_getPreviousDecl(x))
end

function getMostRecentDecl(x::AbstractRecordDecl)
    @check_ptrs x
    return RecordDecl(clang_RecordDecl_getMostRecentDecl(x))
end

function hasFlexibleArrayMember(x::AbstractRecordDecl)
    @check_ptrs x
    return clang_RecordDecl_hasFlexibleArrayMember(x)
end

function isAnonymousStructOrUnion(x::AbstractRecordDecl)
    @check_ptrs x
    return clang_RecordDecl_isAnonymousStructOrUnion(x)
end

function isInjectedClassName(x::AbstractRecordDecl)
    @check_ptrs x
    return clang_RecordDecl_isInjectedClassName(x)
end

function isLambda(x::AbstractRecordDecl)
    @check_ptrs x
    return clang_RecordDecl_isLambda(x)
end

function isCapturedRecord(x::AbstractRecordDecl)
    @check_ptrs x
    return clang_RecordDecl_isCapturedRecord(x)
end

function getDefinition(x::AbstractRecordDecl)
    @check_ptrs x
    return RecordDecl(clang_RecordDecl_getDefinition(x))
end

function isOrContainsUnion(x::AbstractRecordDecl)
    @check_ptrs x
    return clang_RecordDecl_isOrContainsUnion(x)
end


text:
"""
brute force search without skips>1
"""
function naivesearch end

function naivesearch(t::SearchSequence)
    return (t,)
end

function naivesearch(s::SearchSequence, t::SearchSequence, i::Integer)
    naivesearch(s,naivesearch(t),i)
end

function naivesearch(s::SearchSequence,p::Tuple, i::Integer,sv::MaybeVector=nothing)
    (t,) = p
    DOSTATS = !(sv isa Nothing)
    n = sizeof(t)
    m = sizeof(s)

    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return something(findnext(isequal(_nthbyte(t,1)), s, i), 0)
    end

    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

    tlast = _nthbyte(t,n)
    if DOSTATS loops = 0 end
    if DOSTATS bloomtests = 0 end
    if DOSTATS bloomskips = 0 end
    i -= 1
    while i <= w
        if DOSTATS loops += 1 end
        if _nthbyte(s,i+n) == tlast
            # check candidate
            if DOSTATS bloomtests += 1 end # record last byte matches
            j = 1
            while j < n
                if _nthbyte(s,i+j) != _nthbyte(t,j)
                    break
                end
                j += 1
                # match found?
                if j == n
                    if DOSTATS recordcase(sv, loops, 0,0,0,0) end
                    return i+1
                end
            end
        end
        i += 1
    end
    if DOSTATS recordcase(sv, loops, bloomtests,0,0,0) end
    0
end


text:
export readconll, rawtext

function readconll(path, columns=Int[])
    doc = []
    sent = []
    lines = open(readlines, path)
    for line in lines
        line = chomp(line)
        if length(line) == 0
            length(sent) > 0 && push!(doc, sent)
            sent = []
        else
            items = split(line, '\t')
            length(items) > 0 && (items = items[columns])
            push!(sent, items)
        end
    end
    length(sent) > 0 && push!(doc, sent)
    doc
end

function rawtext(filename)
    strs = String[]
    lines = open(readlines, filename)
    for line in lines
        line = chomp(line)
        length(line) == 0 && continue
        items = split(line, '\t')
        tag = items[3]
        if tag != "_"
            for c in tag
                if c == 'S'
                    push!(strs, " ")
                elseif c == 'N'
                    push!(strs, "\n")
                else
                    throw("Invalid char: $(c)")
                end
            end

        end
        push!(strs, items[1])
    end
    join(strs, "")
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
# ====================================================
# Homogeneization of linear systems
# ====================================================

# no-op
homogenize(ivp::IVP{LCS{N, MT}, ST}) where {N, MT<:AbstractMatrix{N}, ST} = ivp
homogenize(sys::LCS{N, MT}) where {N, MT<:AbstractMatrix{N}} = sys

"""
    homogenize(ivp::IVP{CLCCS{N,MT,IdentityMultiple{N},XT,ConstantInput{SI}},ST}) where {N, MT<:AbstractMatrix{N}, XT<:LazySet{N}, SI<:Singleton{N}, ST<:LazySet{N}}

Transform an inhomogeneous linear initial-value problem into an homogeneous one
by introducing auxiliary state variables.

### Input

- `ivp` -- initial-value problem

### Output

Homogeneous initial-value problem.

### Notes

This function transforms the canonical initial-value problem ``x' = Ax + u``,
``x ∈ X`` with ``u(0) ∈ U = {u}`` (a singleton) into an homogeneous problem
without inputs ``y' = Â * y``, ``y ∈ Y``.
"""
function homogenize(ivp::IVP{CLCCS{N,MT,IdentityMultiple{N},XT,ConstantInput{SI}},ST}) where {N, MT<:AbstractMatrix{N}, XT<:LazySet{N}, SI<:Singleton{N}, ST<:LazySet{N}}
    # homogenized state matrix
    U = inputset(ivp) |> ReachabilityAnalysis.next_set
    A = state_matrix(ivp)
    Â = _homogenize_state_matrix(A, U)

    # homogenized input set
    X0 = initial_state(ivp)
    Y0 = _homogenize_initial_state(X0)

    X = stateset(ivp)
    Y = _homogenize_stateset(X)
    ivph = IVP(CLCS(Â, Y), Y0)
    return ivph
end

"""
    homogenize(sys::SOACS)

Transform an inhomogeneous second order system into an homogeneous one
by introducing auxiliary state variables.

### Input

- `sys` -- second order system

### Output

First-order homogeneous system.

### Notes

This function transforms the second-order system ``Mx'' + Cx' + Kx = b`` into a
first-order, homogeneous one, ``y' = Â * y``. It is assumed that the matrix ``M``
is invertible.
"""
function homogenize(sys::SOACS)
    @unpack M, C, K, b = sys

    n = size(M, 1)
    invM = inv(M)
    Zn = spzeros(n, n)
    In = Matrix(1.0I, n, n)

    A = [Zn           In     ;
        -invM*K       -invM*C]
    b0 = vcat(zeros(n), invM * b)

    m = 2n
    Aext = spzeros(m+1, m+1)
    Aext[1:m, 1:m] .= A
    Aext[1:m, m+1] .= b0

    return @system(x' = Aext*x)
end

function _homogenize_state_matrix(A::AbstractMatrix, U::Singleton)
    u = element(U)
    n = size(A, 1)
    Â = zeros(n+1, n+1)
    Â[1:n, 1:n] .= A
    Â[1:n, n+1] .= u
    return Â
end

function _homogenize_initial_state(X0::Singleton{N}) where {N}
    x0 = element(X0)
    y0 = vcat(x0, one(N))
    Y0 = Singleton(y0)
    return Y0
end

function _homogenize_initial_state(X0::LazySet{N}, dims=1) where {N}
    Y0 = X0 × Singleton(ones(N, dims))
    return Y0
end

function _homogenize_stateset(X::Universe, ndims=1)
    return Universe(dim(X) + ndims)
end

function _homogenize_stateset(X::LazySet, ndims=1)
    return X × Universe(ndims)
end

function _homogenize_stateset(X::HalfSpace{N}, ndims=1) where {N}
    a′ = vcat(X.a, fill(zero(N), ndims))
    return HalfSpace(a′, X.b)
end


text:
using Mise
using Test

using DataFrames
using StatsBase: mean, std, mean_and_std, zscore
using Dates, TimeZones
using JuliaDB
using Flux
using CuArrays

CuArrays.allowscalar(false)

include("data.jl")
include("utils.jl")
include("evaluation.jl")


text:
@testset "double pendulum" begin
    Random.seed!(6)
    lc1 = -0.5
    l1 = -1.
    m1 = 1.
    I1 = 0.333 # about joint instead of CoM in URDF
    lc2 = -1.
    l2 = -2.
    m2 = 1.
    I2 = 1.33 # about joint instead of CoM in URDF
    g = -9.81

    axis = SVector(0., 1., 0.)

    double_pendulum = Mechanism(RigidBody{Float64}("world"); gravity = SVector(0, 0, g))
    world = root_body(double_pendulum)

    # create first body and attach it to the world via a revolute joint
    inertia1 = SpatialInertia(CartesianFrame3D("upper_link"), moment=I1 * axis * axis', com=SVector(0, 0, lc1), mass=m1)
    body1 = RigidBody(inertia1)
    joint1 = Joint("shoulder", Revolute(axis))
    joint1_to_world = one(Transform3D, joint1.frame_before, default_frame(world))
    attach!(double_pendulum, world, body1, joint1, joint_pose = joint1_to_world)

    inertia2 = SpatialInertia(CartesianFrame3D("lower_link"), moment=I2 * axis * axis', com=SVector(0, 0, lc2), mass=m2)
    body2 = RigidBody(inertia2)
    joint2 = Joint("elbow", Revolute(axis))
    joint2_to_body1 = Transform3D(joint2.frame_before, default_frame(body1), SVector(0, 0, l1))
    attach!(double_pendulum, body1, body2, joint2, joint_pose = joint2_to_body1)

    @test findbody(double_pendulum, string(body1)) == body1
    @test_throws ErrorException findbody(double_pendulum, "bla")
    @test findjoint(double_pendulum, string(joint2)) == joint2
    @test_throws ErrorException findjoint(double_pendulum, "bla")

    x = MechanismState(double_pendulum)
    rand!(x)

    # from http://underactuated.csail.mit.edu/underactuated.html?chapter=3
    q1 = configuration(x, joint1)[1]
    q2 = configuration(x, joint2)[1]
    v1 = velocity(x, joint1)[1]
    v2 = velocity(x, joint2)[1]
    c1 = cos(q1)
    c2 = cos(q2)
    s1 = sin(q1)
    s2 = sin(q2)
    s12 = sin(q1 + q2)
    T1 = 1/2 * I1 * v1^2
    T2 = 1/2 * (m2 * l1^2 + I2 + 2 * m2 * l1 * lc2 * c2) * v1^2 + 1/2 * I2 * v2^2 + (I2 + m2 * l1 * lc2 * c2) * v1 * v2

    M11 = I1 + I2 + m2 * l1^2 + 2 * m2 * l1 * lc2 * c2
    M12 = I2 + m2 * l1 * lc2 * c2
    M22 = I2
    M = [M11 M12; M12 M22]

    C11 = -2 * m2 * l1 * lc2 * s2 * v2
    C12 = -m2 * l1 * lc2 * s2 * v2
    C21 = m2 * l1 * lc2 * s2 * v1
    C22 = 0
    C = [C11 C12; C21 C22]

    G = [m1 * g * lc1 * s1 + m2 * g * (l1 * s1 + lc2 * s12); m2 * g * lc2 * s12]

    v̇ = similar(velocity(x))
    rand!(v̇)
    τ = inverse_dynamics(x, v̇)
    v = velocity(x)

    @test isapprox(T1, kinetic_energy(x, body1), atol = 1e-12)
    @test isapprox(T2, kinetic_energy(x, body2), atol = 1e-12)
    @test isapprox(M, mass_matrix(x), atol = 1e-12)
    @test isapprox(τ, M * v̇ + C * v + G, atol = 1e-12)

    # compare against URDF
    for revolute_joint_type in [Revolute{Float64}, SinCosRevolute{Float64}]
        double_pendulum_urdf = parse_urdf(joinpath(@__DIR__, "urdf", "Acrobot.urdf"),
            remove_fixed_tree_joints=false, revolute_joint_type=revolute_joint_type)
        x_urdf = MechanismState(double_pendulum_urdf)
        for (i, j) in enumerate(joints(double_pendulum))
            urdf_joints = collect(joints(double_pendulum_urdf))
            index = findfirst(joint -> string(joint) == string(j), urdf_joints)
            j_urdf = urdf_joints[index]
            set_configuration!(x_urdf, j_urdf, configuration(x, j)[1])
            set_velocity!(x_urdf, j_urdf, velocity(x, j))
        end
        v̇ = similar(velocity(x_urdf))
        rand!(v̇)
        τ = inverse_dynamics(x_urdf, v̇)
        urdf_bodies = collect(bodies(double_pendulum_urdf))
        urdf_upper_link = urdf_bodies[findfirst(b -> string(b) == string(body1), urdf_bodies)]
        urdf_lower_link = urdf_bodies[findfirst(b -> string(b) == string(body2), urdf_bodies)]
        @test isapprox(T1, kinetic_energy(x_urdf, urdf_upper_link), atol = 1e-12)
        @test isapprox(T2, kinetic_energy(x_urdf, urdf_lower_link), atol = 1e-12)
        @test isapprox(M, mass_matrix(x_urdf), atol = 1e-12)
        @test isapprox(τ, M * v̇ + C * v + G, atol = 1e-12)
    end
end


text:
function PropagatorsParallelizer(arg0::List, arg1::MultiSatStepHandler)
    return PropagatorsParallelizer((List, MultiSatStepHandler), arg0, arg1)
end

function get_propagators(obj::PropagatorsParallelizer)
    return jcall(obj, "getPropagators", List, ())
end

function propagate(obj::PropagatorsParallelizer, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "propagate", List, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
import Future.randjump

randjump_(rng::MersenneTwister, num) = randjump(rng, num);
randjump_(rng, num) = rng;

inner_seed!(tmp_rng, seed) = Random.seed!(tmp_rng, seed)
#inner_seed!(tmp_rng::VectorizedRNG.AbstractPCG,seed)= return nothing;

function set_seed(mcConfig::MonteCarloConfiguration{type1, type2, type3, type4, type5}) where {type1 <: Integer, type2 <: Integer, type3 <: AbstractMonteCarloMethod, type4 <: SerialMode, type5 <: Random.AbstractRNG}
    tmp_rng = deepcopy(mcConfig.rng)
    inner_seed!(tmp_rng, mcConfig.seed)
    mcConfig.rng = randjump_(tmp_rng, mcConfig.offset)
    Random.seed!(mcConfig.seed * 10 + 10)

    return
end


text:
export ProblemSave, save, loadpb
using JLD
struct ProblemSave

    param_problem :: Dict
    model   :: Symbol
    param_model   :: Dict
    solver  :: Symbol
    param_solver   :: Dict
    data    :: Data

    function ProblemSave(param_p, model, param_m, solver, param_s, data)
         new( param_p, model, param_m, solver, param_s, deepcopy(data))
    end

end

import Base.convert

function convert(::Type{ProblemSave}, p :: Problem )

    model   = Symbol(typeof(p.model))
    solver  = Symbol(typeof(p.solver))

    param_problem = p.param
    if :param in fieldnames(typeof(p.model))
        param_model=p.model.param
    else
        param_model = NamedTuple()
    end
    if :kwargs in fieldnames(typeof(p.model))
        kwargs=zip(keys(p.model.kwargs),values(p.model.kwargs))
        param_model=merge(param_model,(kwargs=kwargs,))
    end
    if :param in fieldnames(typeof(p.solver))
        param_solver=p.solver.param
    else
        param_solver = NamedTuple()
    end
    param_problem = Dict(pairs(param_problem))
    param_model   = Dict(pairs(param_model))
    param_solver  = Dict(pairs(param_solver))
    data    = p.data

    ProblemSave(param_problem, model, param_model, solver, param_solver, data)

end

function convert(::Type{Problem}, p :: ProblemSave)

    param_p = (;p.param_problem...)
    param_m = (;p.param_model...)
    param_s = (;p.param_solver...)

    # Reconstructs the model
    model = try # try with set of parameters defined by the model
        getfield(WaterWaves1D, p.model)(param_m)
    catch m
        if isa(m,AbstractModel) == false # if error, then use set of parameters defined by the problem
            model = getfield(WaterWaves1D, p.model)(param_p)
        end
    end
    if :kwargs in fieldnames(typeof(model))
        model = getfield(WaterWaves1D, p.model)(param_m;param_m.kwargs...)
    end

    # Reconstructs the initial data
    mesh=Mesh(param_p)
    U=first(p.data.U)
    initial = Init(mesh,U[:,1],U[:,2])

    # Reconstructs the solver (may not work with other user-defined solvers)
    solver = try
        getfield(WaterWaves1D, p.solver)(param_s)
    catch s
        if isa(s,TimeSolver) == false
            solver = try
                getfield(WaterWaves1D, p.solver)(model)
            catch s
                if isa(s,TimeSolver) == false
                    solver = try
                        getfield(WaterWaves1D, p.solver)()
                    catch s
                        if isa(s,TimeSolver) == false
                            @warn "Cannot set up the solver. Choose RK4."
                            solver = RK4(model)
                        end
                    end
                end
            end
        end
    end

    # Reconstructs the structure of the problem
    pb = Problem(model, initial, param_p; solver = solver)

    # Writes the raw data
    pb.data = p.data

    return pb

end

import JLD.save

"""
    `save(p::Problem,name::String)`

Saves the content of an object `Problem` into the file `name.jld`.
"""

function save(p::Problem,name::String)

    filename = string(name,".jld")
    JLD.save(filename, name, convert(ProblemSave,p))

end

import JLD.load

export loadpb

"""
    `loadpb(name::String)`

Loads the contents of the file `name.jld` as a problem of type `:Problem`.
"""
function loadpb(name::String)

    convert(Problem, JLD.load(string(name,".jld"), name ))

end


text:
const should_not_warn = Set{Function}()

"""
    @should_not_warn function foo(x)
      ...
    end

Add `foo` to the list of functions in which no warnings may occur (checkd by `@check`).
"""
macro should_not_warn(expr)
  quote
    fun = $(esc(expr))
    push!(should_not_warn, fun)
    fun
  end
end

"""
    check(f::Function)

Run Traceur on `f`, and throw an error if any warnings occur inside functions
tagged with `@should_not_warn`.
"""
function check(f; nowarn=Any[], kwargs...)
  failed = false
  wp = warning_printer()
  result = trace(f; kwargs...) do warning
    ix = findfirst(warning.stack) do call
      call.f in should_not_warn || call.f in nowarn
    end
    if ix != nothing
      tagged_function = warning.stack[ix].f
      message = "$(warning.message) (called from $(tagged_function))"
      warning = Warning(warning.call, warning.line, message, warning.stack)
      wp(warning)
      failed = true
    end
  end
  @assert !failed "One or more warnings occured inside functions tagged with `@should_not_warn`"
  result
end

"""
    @check fun(args...) nowarn=[] maxdepth=typemax(Int)

Run Traceur on `fun`, and throw an error if any warnings occur inside functions
tagged with `@should_not_warn` or specified in `nowarn`.
"""
macro check(expr, args...)
  quote
      check(() -> $(esc(expr)); $(map(esc, args)...))
    end
end


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
theme(:wong)
gr(aspect_ratio=1,legend=:bottomright,axis=nothing,border=nothing,ticks=nothing,legendfontsize=12,xlims=[-3,7],
   markerstrokewidth=0,markersize=1)

trans_x = 1.2
trans_y = 1.6

plot([-2;-2].+2*4/8,[-2;2],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2;-2+trans_x].+2*4/8,[2;2+trans_y],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2+trans_x;-2+trans_x].+2*4/8,[2+trans_y;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x].+2*4/8,[-2;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)

plot!([-2;-2].+3*4/8,[-2;2],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2;-2+trans_x].+3*4/8,[2;2+trans_y],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2+trans_x;-2+trans_x].+3*4/8,[2+trans_y;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x].+3*4/8,[-2;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)

plot!([-2;-2].+5*4/8,[-2;2],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2;-2+trans_x].+5*4/8,[2;2+trans_y],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2+trans_x;-2+trans_x].+5*4/8,[2+trans_y;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x].+5*4/8,[-2;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)

plot!([-2;-2].+6*4/8,[-2;2],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2;-2+trans_x].+6*4/8,[2;2+trans_y],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2+trans_x;-2+trans_x].+6*4/8,[2+trans_y;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x].+6*4/8,[-2;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)

plot!([-2;-2].+7*4/8,[-2;2],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2;-2+trans_x].+7*4/8,[2;2+trans_y],seriescolor=:lightgrey,lw=3,label = nothing)
plot!([-2+trans_x;-2+trans_x].+7*4/8,[2+trans_y;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x].+7*4/8,[-2;-2+trans_y],seriescolor=:lightgrey,lw=3,ls=:dash,label = nothing)

plot!([-2;-2],[-2;2],seriescolor=:dimgrey,lw=3,label = nothing)
plot!([-2;2],[2;2],seriescolor=:dimgrey,lw=3,label = nothing)
plot!([2;-2],[-2;-2],seriescolor=:dimgrey,lw=3,label = nothing)

plot!([-2;-2].+trans_x,[-2;2].+trans_y,seriescolor=:dimgrey,lw=3,ls=:dash,label = nothing)
plot!([-2;2].+trans_x,[2;2].+trans_y,seriescolor=:dimgrey,lw=3,label = nothing)
plot!([2;-2].+trans_x,[-2;-2].+trans_y,seriescolor=:dimgrey,lw=3,ls=:dash,label = nothing)

plot!([-2;-2+trans_x],[-2;-2+trans_y],seriescolor=:dimgrey,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x],[2;2+trans_y],seriescolor=:dimgrey,lw=3,label = nothing)

# color1 = :royalblue
# color2 = :deepskyblue
# color3 = :lightskyblue1

color1 = :firebrick4
color2 = :firebrick2
color3 = :coral

# First Fourier slice
plot!([-2;-2].+4/8,[-2;2],seriescolor=color1,lw=3,label = nothing)
plot!([-2;-2+trans_x].+4/8,[2;2+trans_y],seriescolor=color1,lw=3,label = "Fourier mode 1")
plot!([-2+trans_x;-2+trans_x].+4/8,[2+trans_y;-2+trans_y],seriescolor=color1,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x].+4/8,[-2;-2+trans_y],seriescolor=color1,lw=3,ls=:dash,label = nothing)

# Second Fourier slice
plot!([-2;-2].+4*4/8,[-2;2],seriescolor=color2,lw=3,label = nothing)
plot!([-2;-2+trans_x].+4*4/8,[2;2+trans_y],seriescolor=color2,lw=3,label = "Fourier mode 4")
plot!([-2+trans_x;-2+trans_x].+4*4/8,[2+trans_y;-2+trans_y],seriescolor=color2,lw=3,ls=:dash,label = nothing)
plot!([-2;-2+trans_x].+4*4/8,[-2;-2+trans_y],seriescolor=color2,lw=3,ls=:dash,label = nothing)

# Third Fourier slice
plot!([2;2+trans_x],[-2;-2+trans_y],seriescolor=color3,lw=3,label = nothing)
plot!([2;2+trans_x],[2;2+trans_y],seriescolor=color3,lw=3,label = nothing)
plot!([2;2],[2;-2],seriescolor=color3,lw=3,label = nothing)
plot!([2;2].+trans_x,[2;-2].+trans_y,seriescolor=color3,lw=3,label = "Fourier mode 8")

plot!([-1.5,1.5],[-2.2,-2.2],arrow=0.4,seriescolor=:black,lw=2.5,label=nothing)
annotate!(0,-2.5,"z direction")
savefig("domain.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
fn = normpath(joinpath(dirname(@__FILE__),"center_sizes.jld2"))
d = JLD2.jldopen(fn,"r")

tmp = d["ctrs"]
ctrs = [SVector(q...) for q in tmp]
rads = d["rads"]

tree = CD.Octree(ctrs, rads)

# extract all the triangles that (potentially) intersect octant (+,+,+)
pred(i) = all(ctrs[i].+rads[i] .> 0)
bb = SVector(0.5, 0.5, 0.5), 0.5
ids = collect(CD.searchtree(pred, tree, bb))
@test length(ids) == 178

N = 100
using DelimitedFiles
buf = readdlm(joinpath(@__DIR__,"assets","ctrs.csv"))
ctrs = vec(collect(reinterpret(SVector{3,Float64}, buf')))
rads = vec(readdlm(joinpath(@__DIR__,"assets","rads.csv")))

tree = CD.Octree(ctrs, rads)
pred(i) = all(ctrs[i] .+ rads[i] .> 0)
bb = @SVector[0.5, 0.5, 0.5], 0.5
ids = collect(CD.searchtree(pred, tree, bb))

@show ids
ids2 = findall(i -> all(ctrs[i].+rads[i] .> 0), 1:N)
@test length(ids2) == length(ids)
@test sort(ids2) == sort(ids)
@test ids == [26, 46, 54, 93, 34, 94, 75, 23, 86, 57, 44, 40, 67, 73, 77, 80]

max_stars_repo_name: krcools/CollisionDetection

max_issues_repo_name: krcools/CollisionDetection

max_forks_repo_name: krcools/CollisionDetection


text:
function get_df(fp::FairnessProblem, models_fn)

	Logging.disable_logging(LogLevel(0))

	#aliasing
	X = fp.task.X
	y = fp.task.y
	measures=fp.measures
	protected_attr = fp.task.grp
	debiasmeasures = fp.task.debiasmeasures
	refGrp = fp.refGrp
	runs=fp.repls
	nfolds=fp.nfolds
	random_seed = fp.seed
	models, model_names = models_fn(fp)

	Logging.disable_logging(LogLevel(0))
	# Here runs implies number of runs of n_folds cv-folds
	df = DataFrame([String[], Float64[],
			[Float64[] for i in 1:length(measures)]..., String[], Int[], Int[]],
			["model", "accuracy",
			[string(measure)*"_disparity" for measure in measures]...,
			"domain", "fold", "replication"]) #domain tells whether result is for in-sample or out-sample
	Random.seed!(random_seed)
	seeds = abs.(rand(Int, runs))
	if model_names==nothing model_names = string.(models) end
	progress = Progress(runs*length(models)*nfolds, 1)
	for i in 1:runs
		cv = StratifiedCV(nfolds=nfolds, shuffle=true, rng=seeds[i])
		tr_tt_pairs = MLJBase.train_test_pairs(cv, 1:length(y),
			categorical(string.(X[!, protected_attr]) .* "-" .* string(y)))
		for i_model in 1:length(models)
			model = models[i_model]
			for j in 1:nfolds
				Random.seed!(seeds[i])
				mach = machine(model, X, y)
				@suppress fit!(mach, rows=tr_tt_pairs[j][1], verbosity=0)

				# Add out-sample performance measures
				ŷ = MMI.predict(mach, rows=tr_tt_pairs[j][2])
				if typeof(ŷ[1])<:MLJBase.UnivariateFinite ŷ = StatsBase.mode.(ŷ) end
				ft = fair_tensor(ŷ, y[tr_tt_pairs[j][2]],
										X[tr_tt_pairs[j][2], protected_attr])
				accVal = accuracy(ft)
				# disps = disparity(measures, ft, refGrp=refGrp)
				# For now, the value of Disparity I will consider will be:
				# (Overall Fairness Value)/(Fairness value for reference Group)
				push!(df, [model_names[i_model], accVal,
				[measure(ft, grp="0")/measure(ft, grp="1") for measure in measures]...,
				"test", j, i])

				# Add in-sample performance measures
				ŷ = MMI.predict(mach, rows=tr_tt_pairs[j][1])
				if typeof(ŷ[1])<:MLJBase.UnivariateFinite ŷ = StatsBase.mode.(ŷ) end
				ft = fair_tensor(ŷ, y[tr_tt_pairs[j][1]],
										X[tr_tt_pairs[j][1], protected_attr])
				accVal = accuracy(ft)
				# disps = disparity(measures, ft, refGrp=refGrp)
				# For now, the value of Disparity I will consider will be:
				# (Overall Fairness Value)/(Fairness value for reference Group)
				push!(df, [model_names[i_model], accVal,
				[measure(ft, grp="0")/measure(ft, grp="1") for measure in measures]...,
				"train", j, i])
				next!(progress)
			end
		end
	end
	CSV.write(joinpath(pwd(), fp.name*"-results.csv"), df)
	return df
end

function get_pareto_df(fp::FairnessProblem, models_fn, alphas=0:0.1:1)

	Logging.disable_logging(LogLevel(0))

	#aliasing
	X = fp.task.X
	y = fp.task.y
	measures=fp.measures
	protected_attr = fp.task.grp
	debiasmeasures = fp.task.debiasmeasures
	refGrp = fp.refGrp
	runs=fp.repls
	nfolds=fp.nfolds
	random_seed = fp.seed
	models, model_names = models_fn(fp, 0.5)

	# Here runs implies number of runs of n_folds cv-folds
	df = DataFrame([String[], Float64[], Float64[],
			[Float64[] for i in 1:length(measures)]..., String[], Int[], Int[]],
			["model", "alpha", "accuracy",
			[string(measure)*"_disparity" for measure in measures]...,
			"domain",  "fold", "replication"])
			#domain tells whether result is for in-sample or out-sample
	Random.seed!(random_seed)
	seeds = abs.(rand(Int, runs))
	if model_names==nothing model_names = string.(models) end
	progress = Progress(length(alphas)*runs*length(models)*nfolds, 1)
	for alpha in alphas
		models, _ = models_fn(fp, alpha)
		for i in 1:runs
			cv = StratifiedCV(nfolds=nfolds, shuffle=true, rng=seeds[i])
			tr_tt_pairs = MLJBase.train_test_pairs(cv, 1:length(y),
				categorical(string.(X[!, protected_attr]) .* "-" .* string(y)))
			for i_model in 1:length(models)
				model = models[i_model]
				for j in 1:nfolds
					Random.seed!(seeds[i])
					mach = machine(model, X, y)
					@suppress fit!(mach, rows=tr_tt_pairs[j][1])
					# Add out-sample performance measures
					ŷ = MMI.predict(mach, rows=tr_tt_pairs[j][2])
					if typeof(ŷ[1])<:MLJBase.UnivariateFinite ŷ = StatsBase.mode.(ŷ) end
					ft = fair_tensor(ŷ, y[tr_tt_pairs[j][2]],
											X[tr_tt_pairs[j][2], protected_attr])
					accVal = accuracy(ft)
					# disps = disparity(measures, ft, refGrp=refGrp)
					# For now, the value of Disparity I will consider will be:
					# (Overall Fairness Value)/(Fairness value for reference Group)
					push!(df, [model_names[i_model], alpha, accVal,
					[measure(ft, grp="0")/measure(ft, grp="1") for measure in measures]...,
					"test", j, i])

					# Add in-sample performance measures
					ŷ = MMI.predict(mach, rows=tr_tt_pairs[j][1])
					if typeof(ŷ[1])<:MLJBase.UnivariateFinite ŷ = StatsBase.mode.(ŷ) end
					ft = fair_tensor(ŷ, y[tr_tt_pairs[j][1]],
											X[tr_tt_pairs[j][1], protected_attr])
					accVal = accuracy(ft)
					# disps = disparity(measures, ft, refGrp=refGrp)
					# For now, the value of Disparity I will consider will be:
					# (Overall Fairness Value)/(Fairness value for reference Group)
					push!(df, [model_names[i_model], alpha, accVal,
					[measure(ft, grp="0")/measure(ft, grp="1") for measure in measures]...,
					"train", j, i])
					next!(progress)
				end
			end
		end
	end
	CSV.write(joinpath(pwd(), fp.name*"-pareto-results.csv"), df)
	return df
end

# Core idea behind this function is to get the output after evaluation, hypothesis testing, etc. that can be
# directly passed to plotting functions or any benchmark related function without need to pass any other argument.
"""
    fairevaluate(classifiers, X, y; measures=nothing, measure=nothing, grp=:class, priv_grps, random_seed=12345, n_grps=6)

Performed paired t-test for each pair of classifier in classifiers and return p values and t statistics.
# Arguments
- `classifiers`: Array of classifiers to compare
- `X`: DataFrame with features and protected attribute
- `y`: Binary Target Variable
- `measures=nothing`: The measures to be evaluated and used for HypothesisTests.
                If this is not specified, the `measure` argument is used
- `measure=nothing`: The performance/fairness measure used to perform hypothesis tests.
                    If no values for measure is passed, then Disparate Impact will be used by default.
- `grp=:class`: Protected Attribute Name
- `priv_grps=nothing`: If default measure i.e. Disparate Impact is used, then pass an array of groups which are privileged in dataset.
- `random_seed=12345`: Random seed to ensure reproducibility
- `n_grps=6`: Number of folds for cross validation
# Returns
A dictionary with following keys vs values is returned
- `measures`: names of the measures
- `classifier_names`: names of the classifiers. If a pipeline is used, it will show pipeline and associated number.
- `results`: 3-dimensional array with evaluation result. Its size is measures x classifiers x fold_number.
- `pvalues`: 3-dimensional array with pvalues for each pair of classifier. Its size is measures x classifiers x classifiers.
- `tstats`:3-dimensional array with tstats for each pair of classifier. Its size is measures x classifiers x classifiers.
"""
function fairevaluate(
    classifiers::Array{<:MLJBase.Model,1}, X, y;
    measure = nothing,
    measures = nothing,
    grp = :class,
    priv_grps = nothing,
    random_seed::Int = 12345,
    n_folds = 6,
    classifier_names = nothing
)
    Random.seed!(random_seed)
    y = coerce(y, OrderedFactor)

    @assert(!(measure==nothing && measures==nothing && priv_grps==nothing))
    if priv_grps!=nothing
        @assert(all(
            priv_grp in levels(X[!, grp]) for priv_grp in priv_grps
            for classifier in classifiers
        ))
        measure = measure == nothing ? DisparateImpact(grp, priv_grps) : measure
    end
    @assert(all(
        target_scitype(classifier) <: AbstractVector{<:Finite}
        for classifier in classifiers
    ))

    n = length(classifiers)

    if measures==nothing measures=[measure] end
    n_measures = length(measures)

    for i in 1:n_measures
        if typeof(measures[i]) <: MetricWrapper
            measures[i].grp = grp
        end
    end

    results = zeros(n_measures, n, n_folds)
    cv = CV(nfolds = n_folds, shuffle=false, rng=random_seed)
    for i = 1:n
        Random.seed!(random_seed)
        operation = istype(classifiers[i], Probabilistic) ? MLJBase.predict_mode : MLJBase.predict
        result = evaluate(
            classifiers[i],
            X,
            y,
            resampling = cv,
            measures = measures,
            operation = operation,
            verbosity=0,
        )
        for j in 1:n_measures
            results[j, i, :] = result.per_fold[j]
        end
    end
    pvalues, tstats = zeros(n_measures, n, n), zeros(n_measures, n, n)
    for k in 1:n_measures
        for i = 1:n
            for j = 1:n
                ttestResult = OneSampleTTest(results[k, i, :], results[k, j, :])
                tstats[k, i, j] = ttestResult.t
                pvalues[k, i, j] = pvalue(ttestResult)
            end
        end
    end
    if classifier_names == nothing
        classifier_names = string.(classifiers)
    end
    dict = Dict()
    dict["measures"] = string.(measures)
    dict["classifier_names"] = classifier_names
    dict["results"] = results
    dict["pvalues"] = pvalues
    dict["tstats"] = tstats
    return dict
end


text:
abstract type stType end
struct fullType <: stType end
struct decreasingType <: stType end

(::Colon)(a::Tuple, b::Tuple) = (:).(a, b)
(::Colon)(a::Tuple, b::Tuple, c::Tuple) = (:).(a, b, c)

@doc """
      n, q, dataSizes, outputSizes, resultingSize = calculateSizes(layers::stParallel{K,1}, outputSubsample, Xsize; totalScales = [-1 for i=1:depth(layers)+1]) where {K}

  * n is a list of the sizes of a single example. It is (m+1)×dataDim
one dimension at a time
  * q is a list of the number of scales, as used in 1D
  * datasizes is a list of tuples of sizes encountered in each layer
  * outputsizes is the same for the output arrays
  * resultingsize is a list of the output size as a result of subsampling
"""
function calculateSizes(layers::stParallel{K,1}, outputSubsample, Xsize; totalScales = [-1 for i = 1:depth(layers)+1], percentage = 0.9) where {K}
    n = Int.(sizes(bsplineType(), layers.subsampling, Xsize[1:1]))
    q = getQ(layers, n, totalScales)
    dataSizes = [Int.([n[i], q[i], Xsize[2:end]...]) for i = 1:depth(layers)+1]
    outputSizes = [Int.([n[i+1], q[i], Xsize[2:end]...]) for i = 1:depth(layers)+1]
    if outputSubsample[1] > 1
        resultingSize = zeros(Int, depth(layers) + 1)
        resultingSubsampling = zeros(depth(layers) + 1)
        # subsampling limited by the second entry
        for (i, x) in enumerate(outputSizes)
            proposedLevel = floor(Int, x[1] / outputSubsample[1])
            if proposedLevel < outputSubsample[2]
                # the result is smaller than the floor
                proposedLevel = outputSubsample[2]
            end
            resultingSize[i] = Int(proposedLevel)
        end
    elseif outputSubsample[2] > 1
        resultingSize = outputSubsample[2] * ones(Int64, depth(layers) + 1)
    else
        resultingSize = [x[1] for x in outputSizes]
    end
    return (n, q, dataSizes, outputSizes, resultingSize)
end

"""

get the sizes of a single example in each layer
"""
function getNs(Xsize, layers::stParallel{K,1}) where {K}
    return Int.(sizes(bsplineType(), layers.subsampling, Xsize[(end):end]))
end

@doc """
  q = getQ(layers, n, totalScales; product=true)
calculate the total number of entries in each layer
"""
function getQ(layers::stParallel{K,1}, n, totalScales; product = true) where {K}
    # first just the number of new scales
    q = [numScales(layers.shears[i], n[i]) - 1 for i = 1:depth(layers)+1]
    # then a product over all previous
    q = [(isnan(totalScales[i]) || totalScales[i] <= 0) ? q[i] : totalScales[i] for i = 1:depth(layers)+1]
    if product
        q = [prod(q[1:i-1]) for i = 1:depth(layers)+1]
        return q
    else
        return [1; reverse(q)[2:end]...]
    end
end

getQ(layers, Xsize; product = true) = getQ(layers, Int.(sizes(bspline, layers.subsampling, Xsize[(end):end])), [-1 for i = 1:depth(layers)+1]; product = product)

"""

a list of the paths in a given transform. E.g. if the first layer has 13 scales and the second layer has 9, then this gives [[], [13], [9,13]].
"""
function getListOfScaleDims(layers, n, totalScales = [-1 for i = 1:depth(layers)+1])
    nScalesEachLayer = getQ(layers, n, totalScales; product = false)
    addedLayers = [nScalesEachLayer[end+2-i:end] for i = 1:depth(layers)+1]
end
"""
return a list of locations for iterating over layer i in the same order as the
transform (depth first, so early layers are the slow changing variable)
"""
function indexOverScales(addedLayers, i)
    if i != 1
        axe = [1:nScales for nScales in addedLayers[i]]
        Iterators.product(axe...)
    else
        return 1
    end
end

function linearFromTuple(js, nScales)
    if length(js) > 1
        js[1] - 1 + nScales[1] * linearFromTuple(js[2:end], nScales[2:end])
    else
        return js[1] - 1
    end
end
function getLastScales(dataDim, daughters, shears)
    size(daughters, 2)
end

"""
    createFFTPlans(layers::stParallel{K, N}, dataSizes; verbose=false, T=Float32, iscomplex::Bool=false)

to be called after distributed. returns a 2D array of futures; the i,jth entry has a future (referenced using fetch) for the fft to be used by worker i in layer j.
"""
function createFFTPlans(layers::stParallel{<:Any,N}, dataSizes;
    T = Float32, iscomplex::Bool = false) where {N}
    @debug "starting to create plans at all"
    @debug "" nwork = nworkers()
    nPlans = getNumPlans(layers)
    FFTs = Array{Future,2}(undef, nworkers(), depth(layers) + 1)
    for i in 1:nworkers(), j in 1:depth(layers)+1
        @debug "i=$(i), j=$(j)"
        if length(layers.n) >= 2
            FFTs[i, j] = remotecall(createRemoteFFTPlan, i,
                dataSizes[j][[1 3:end]], T,
                iscomplex, nPlans)
        else
            fftPlanSize = (dataSizes[j][1] * 2, dataSizes[j][3:end]...)
            # in the last layer
            @debug "plans are going to be size $(fftPlanSize)"
            FFTs[i, j] = remotecall(createRemoteFFTPlan, i, fftPlanSize,
                T, iscomplex, nPlans)
        end
    end

    return FFTs
end

# only in the 1D case where we're using either Morlet or Paul wavelets do we
# need both a rfft and an fft
function getNumPlans(layers::stParallel{<:Any,1})
    if typeof(layers.shears[1].waveType) <: Union{Morlet,Paul}
        return 2
    else
        return 1
    end
end

"""
the 1D version
"""
function createRemoteFFTPlan(sizeOfArray, T, iscomplex, nPlans::Int)
    if iscomplex
        return plan_fft!(zeros(Complex{T}, sizeOfArray...), (1,), flags = FFTW.PATIENT)
    elseif nPlans == 1
        return plan_rfft(zeros(T, sizeOfArray...), (1,), flags = FFTW.PATIENT)
    elseif nPlans == 2
        return (plan_rfft(zeros(T, sizeOfArray...), (1,), flags = FFTW.PATIENT),
            plan_fft!(zeros(Complex{T}, sizeOfArray...), (1,), flags = FFTW.PATIENT))
    end
end

function remoteMult(x::Future, y)
    fetch(x) * y
end
remoteMult(y, x::Future) = remoteMultiply(x, y)

function remoteDiv(x::Future, y)
    fetch(x) / y
end
remoteDiv(y, x::Future) = remoteMultiply(x, y)

function computeAllWavelets(layers, N; outputSubsample = (-1, -1), totalScales = [-1 for i = 1:depth(layers)+1], percentage = 0.9)
    n, q, dataSizes, outputSizes, resultingSize = calculateSizes(layers, outputSubsample, N, totalScales = totalScales, percentage = percentage)
    daughters, ω = computeWavelets(n[1], layers.shears[1]; nScales = -1)
    allDaughters = Array{Array{eltype(daughters),2},1}(undef, length(n))
    allDaughters[1] = daughters
    for (ii, x) in enumerate(n[1:end-1])
        daughters, ω = computeWavelets(n[ii],
            layers.shears[ii];
            nScales = -1)
        allDaughters[ii] = daughters
    end
    daughters, ω = computeWavelets(n[end],
        layers.shears[end];
        nScales = 0)
    allDaughters[end] = daughters
    return allDaughters
end
function plotAllWavelets(allDaughters)
    return plot(plot([allDaughters[1] sum(allDaughters[1], dims = 2)], legend = false, title = "layer 1"),
        heatmap([allDaughters[1] sum(allDaughters[1], dims = 2)]),
        plot([allDaughters[2] sum(allDaughters[2], dims = 2)], legend = false, title = "layer 2"),
        heatmap([allDaughters[2] sum(allDaughters[2], dims = 2)]),
        plot(allDaughters[3][:, 1], legend = false, title = "output last"),
        heatmap(allDaughters[3][:, 1:1]),
        layout = (3, 2)
    )
end


text:
module PasteMyst

import JSON
import HTTP

include("pastestruct.jl")

function getRequester(URI::String)
	r = HTTP.get(URI)
	if r.status != 200
		return r.status
	end
	return JSON.parse(String(r.body))
end

# /paste (GET)

function getPaste(pasteID::String)
	"""
	Get a Paste from the ID specified.
	"""
	return getRequester(string("https://paste.myst.rs/api/v2/paste/", pasteID))
end

# /data (GET)

function getLanguageInfoByExtension(extension::String)
	"""
	Get a Data about the Language from Extension.
	"""
	return getRequester(string("https://paste.myst.rs/api/v2/data/languageExt?extension=", extension))
end

function getLanguageInfoByName(name::String)
	"""
	Get a Data about the Language by Name.
	Note that the Name must be Percent encoded.
	"""
	return getRequester(string("https://paste.myst.rs/api/v2/data/language?name=", name))
end

# /time (GET)

function expiresInToEpoch(createdAt::UInt64, expiresIn::String)
	"""
	Get a Unix Timestamp (Unix Epoch) from the expiresIn value.
	"""
	return getRequester(string("https://paste.myst.rs/api/v2/time/expiresInToUnixTime?createdAt=", createdAt, "&expiresIn=", expiresIn))["result"]
end

# /user (GET)

function isUserExist(username::String)
	"""
	Check If a User exist.
	This method returns boolean.
	Or an Error code If there was an Error sending a Request to PasteMyst.
	Note that the User must be marked as public.
	"""
	r = HTTP.request("GET", string("https://paste.myst.rs/api/v2/user/", username, "/exists"))
	if r.status == 200
		return true
	elseif r.status == 404
		return false
	else
		return r.status
	end
end

function getUserInfo(username::String)
	"""
	Get an Info about the User from Username.
	"""
	return getRequester(string("https://paste.myst.rs/api/v2/user/", username))
end

function getCurrentUserInfo(authkey::String)
	"""
	Get a User info from the Authorization key (API key) provided.
	"""
	r = HTTP.request("GET", "https://paste.myst.rs/api/v2/user/self", ["Authorization"=>authkey])
	if r.status != 200
		return r.status
	end
	return JSON.parse(String(r.body))
end

function getUserPastes(authkey::String)
	"""
	Get the List of all Pastes the user provided has created.
	Note that the User must provide as an API key.
	"""
	r = HTTP.request("GET", "https://paste.myst.rs/api/v2/user/self/pastes", ["Authorization"=>authkey])
	if r.status != 200
		return r.status
	end
	return JSON.parse(String(r.body))
end

# /paste (POST)

function pastyCreateInfoToDict(createInfo::PastyCreateInfo)
    """
    Convert a Single PastyCreateInfo into a Dictionary
    """
	return Dict(
		"title"=>createInfo.title,
		"language"=>createInfo.language,
		"code"=>createInfo.code
	)
end

function pastiesCreateInfoToDict(createInfo::Vector{PastyCreateInfo})
	"""
	Convert all PastyCreateInfo into a Dictionary.
	"""
	pasties::Vector{Dict} = []
	foreach(i -> pasties = append!(pasties, [pastyCreateInfoToDict(i)]), createInfo)
	return pasties
end

function createInfoToDict(createInfo::PasteCreateInfo)
	"""
	Converts the PasteCreateInfo struct info a Dictionary
	"""
	return Dict(
		"title"=>createInfo.title,
		"expiresIn"=>createInfo.expiresIn,
		"isPrivate"=>createInfo.isPrivate,
		"isPublic"=>createInfo.isPublic,
		"tags"=>createInfo.tags,
		"pasties"=>pastiesCreateInfoToDict(createInfo.pasties)
	)
end

function createPaste(createInfo::PastyCreateInfo)
	"""
	This method create a Paste on PasteMyst website without a Authorization key (API key).
	Returns:
		The paste object If the Paste is successfully created.
		Else return false
	"""
	if (createInfo.Tags != "")
		return false
	end

	r = HTTP.request("POST", "https://paste.myst.rs/api/v2/paste/",
		[], string(JSON.json(createInfoToDict(PasteCreateInfo)))
	)
	if r.status != 200
		return false
	end
	return JSON.parse(String(r.body))
end

function dictToPasty(pastyinfo::Dict)
	"""
	Converts a Dictionary to Pasty (struct).
	"""
	return pasty(
		pastyinfo["_id"],
		pastyinfo["language"],
		pastyinfo["title"],
		pastyinfo["code"]
	)
end

function listToPastiesList(pasties::Vector)
	"""
	Convert a List of Pasty (Dictionary) to a List of Pasty (struct).
	"""
	out::Vector{pasty} = []
	foreach(i -> out = append!(out, [dictToPasty(i)]), pasties)
	return out
end

function dictToEdit(editinfo::Dict)
	"""
	Converts a Dictionary to Pasty (struct).
	"""
	return edit(
		editinfo["_id"],
		editinfo["editId"],
		editinfo["editType"],
		editinfo["metadata"],
		editinfo["edit"],
		editinfo["editedAt"]
	)
end

function listToEditList(edits::Vector)
	"""
	Convert a List of Pasty (Dictionary) to a List of Pasty (struct).
	"""
	out::Vector{edit} = []
	foreach(i -> out = append!(out, [dictToEdit(i)]), edits)
	return out
end

function dictToPaste(pasteinfo::Dict)
	"""
	Converts a Dicitionary into a paste struct.
	"""
	return paste(
		pasteinfo["_id"],
		pasteinfo["ownerId"],
		pasteinfo["title"],
		pasteinfo["createdAt"],
		pasteinfo["expiresIn"],
		pasteinfo["deletesAt"],
		pasteinfo["stars"],
		pasteinfo["isPrivate"],
		pasteinfo["isPublic"],
		pasteinfo["encrypted"],
		pasteinfo["tags"],
		listToPastiesList(pasteinfo["pasties"]),
		listToEditList(pasteinfo["edits"]),
	)
end

function createPaste(createInfo::PasteCreateInfo, authkey::String)
	"""
	This method create a Paste on PasteMyst website with an Authorization key (API key).
	Returns:
		The paste object If the Paste is successfully created.
		Else return false
	"""
	r = HTTP.request("POST", "https://paste.myst.rs/api/v2/paste/",
		["Authorization"=>authkey, "content-type"=>"application/json"],
		string(JSON.json(createInfoToDict(createInfo)))
	)
	if r.status != 200
		return false
	end
	return dictToPaste(JSON.parse(String(r.body)))
end

# /paste/ (DELETE)

function deletePaste(authkey::String, pasteID::String)
	"""
	This function delete a Paste on PasteMyst.
	This method requires the Authorization key (API key) to work.
	"""
	r = HTTP.request("DELETE", string("https://paste.myst.rs/api/v2/paste/", pasteID), ["Authorization"=>authkey])
	if r.status != 200
		return r.status
	end
end

# /paste/ (PATCH)

function pastyInfoToDict(pastyinfo::pasty)
	"""
	Converts a Pasty struct into a dictionary.
	"""
	return Dict(
		"_id"=>pastyinfo._id,
		"title"=>pastyinfo.title,
		"language"=>pastyinfo.language,
		"code"=>pastyinfo.code
	)
end

function pastiesInfoToDict(pasties::Vector{pasty})
	"""
	Converts a List of Pasty info a List of Pasties Dictionary.
	"""
	out::Vector{Dict} = []
	foreach(i -> out = append!(out, [pastyInfoToDict(i)]), pasties)
	return out
end

function pasteToDict(pasteinfo::paste)
	"""
	Converts a paste struct into a Dictionary.
	"""
	return Dict(
		"_id"=>pasteinfo._id,
		"ownerId"=>pasteinfo.ownerId,
		"title"=>pasteinfo.title,
		"createdAt"=>pasteinfo.createdAt,
		"expiresIn"=>pasteinfo.expiresIn,
		"deletesAt"=>pasteinfo.deletesAt,
		"stars"=>pasteinfo.stars,
		"isPrivate"=>pasteinfo.isPrivate,
		"isPublic"=>pasteinfo.isPublic,
		"encrypted"=>pasteinfo.encrypted,
		"tags"=>pasteinfo.tags,
		"pasties"=>pastiesInfoToDict(pasteinfo.pasties),
		"edits"=>pasteinfo.edits
	)
end

function pasteEditInfoToDict(pasteinfo::pasteEditInfo)
	"""
	This function converts pasteEditInfo into a Dictionary.
	"""
	return Dict(
		"title"=>pasteinfo.title,
		"isPrivate"=>pasteinfo.isPrivate,
		"isPublic"=>pasteinfo.isPublic,
		"tags"=>pasteinfo.tags,
		"pasties"=>pastiesInfoToDict(pasteinfo.pasties)
	)
end

function pasteToPasteEditInfo(pasteinfo::paste)
	"""
	This function converts paste into pasteEditInfo.
	"""
	return pasteEditInfo(
		pasteinfo.title,
		pasteinfo.isPrivate,
		pasteinfo.isPublic,
		pasteinfo.tags,
		pasteinfo.pasties
	)
end

function editPaste(authkey::String, pasteid::String, pasteinfo::pasteEditInfo)
	"""
	This function edit a Paste on Pastemyst.
	An Authorization key (API key) is required.
	Returns:
		Return the Returned JSON body as Paste struct.
		Else the HTTP error code.
	"""
	r = HTTP.request(
		"PATCH",
		string("https://paste.myst.rs/api/v2/paste/", pasteid),
		["Authorization"=>authkey, "content-type"=>"application/json"],
		string(JSON.json(pasteEditInfoToDict(pasteinfo)))
	)
	if r.status != 200
		return r.status
	end
	return dictToPaste(JSON.parse(String(r.body)))
end

end


text:
"""Splits a global data file into a number of chunks that are later processed in parallel"""
function split_raw_data(j)
  # todo: potentially add caching. only execute if input filename has been modified after the splitdir contents
  for filename in j.inputs
    wrkrs = workers()
    input_filename = joinpath(datadir, filename)
    !isfile(input_filename) && return true # not in data/ but already in split/
    info("Splitting $input_filename using $(length(wrkrs)) workers")
    output_filenames = Array{String}(length(wrkrs))
    for i in 1:length(wrkrs)
      output_dir = joinpath(splitdir, string(wrkrs[i]))
      !isdir(output_dir) && mkdir(output_dir)
      output_filenames[i] = joinpath(output_dir, filename)
    end
    # open the files for writing (1 file per core)
    handles = [open(f, "w") for f in output_filenames]
    handle = 1
    for ln in readlines(input_filename)
      current_output_file = output_filenames[handle]
      # println("current handle: $(handle) and current output filename $(output_filenames[handle])")
      current_handle = handles[handle]
      write(current_handle, ln * "\n") # write line into current handle
      handle += 1 # get next handle
      if handle > length(handles)
        handle = 1 # cycle back to first handle once handles have been exhausted
      end
    end
    map(close, handles) # close all output handles
  end
  return true
end

job_output_dir(j) = job_output_dir(j.jobid)
job_output_dir(jobid::Int) = joinpath(outputdir, string(jobid))

"""Saves the output of a job in its respective output directory"""
function save_job_output(jobid::Int, header::String, output)
  output_file = joinpath(job_output_dir(jobid), "output.csv")
  isfile(output_file) && rm(output_file)
  fh = open(output_file, "w")
  header != "" && write(fh, header * "\n")
  for s in output
    write(fh, s * "\n")
  end
  close(fh)
  true
end

"""Creates an output area for a particular job. If such an area already exists, it is removed and recreated"""
function create_job_area(j)
  info("Creating job area...")
  joboutputdir = job_output_dir(j)
  isdir(joboutputdir) && rm(joboutputdir, force=true, recursive=true)
  mkdir(joboutputdir)
end

"""Reads an inputfile for a particular job"""
function read_sink(j, phase)
  inputfile = joinpath(job_output_dir(j), string(myid()) * "." * phase)
  open(inputfile)
end

"""Writes to an outputfile of a particular job"""
function write_sink(j)
  outputfile = joinpath(splitdir, string(myid()), j.outputfilename)
  isfile(outputfile) && rm(outputfile)
  open(outputfile, "w")
end

function clean_split_dir()
  rm(splitdir, force=true, recursive=true)
  mkdir(splitdir)
end


text:
## In this notebook, we'll look at least squares on biometric dataset
# with a stochastic gradient descent and least-squares angle.
##
using Plots
using ORCA
theme(:dark)
##
# https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv
using CSV
df = CSV.read(download("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv"))
X = Float64.([df[:weight] df[:height]])
# Filter the data
goodpts = X[:,1] .<= 115
Xf = X[goodpts,:]
##
x = Xf[:,1]
A = [x ones(length(x))]
b = Xf[:,2]
ab = A\b # solve the least squares problem
## Show k random samples of the quadratic surface
using Random## In this notebook, we'll look at least squares on biometric dataset
# with a stochastic gradient descent and least-squares angle.
##
using Plots
using ORCA
theme(:dark)
##
# https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv
using CSV
df = CSV.read(download("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv"))
X = Float64.([df[:weight] df[:height]])
# Filter the data
goodpts = X[:,1] .<= 115
Xf = X[goodpts,:]
##
x = Xf[:,1]
A = [x ones(length(x))]
b = Xf[:,2]
ab = A\b # solve the least squares problem
## Show k random samples of the quadratic surface
using Random
using LinearAlgebra
plotlyjs()
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)/size(A,1)
ezsurf!(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot!(x, y, Z, st=:surface, alpha=0.4)
end

## Show the same thing with a log10
quadratic_leastsq(A,b) = x-> log10(0.5*norm(A*x - b)^2)/size(A,1)

Random.seed!(0)
nsamp = 10
k = 25

plot(framestyle=:none,colorbar=false,size=(300,300),dpi=300)
plot!(colorbar=false,size=(900,900))

Random.seed!(0)
for j=1:k
  subset = randperm(size(A,1))[1:nsamp] # choose 25 of 200 samples.
  ezsurf!(-3:0.05:3,-50:0.25:250,
    quadratic_leastsq(A[subset,:],b[subset,:]))
end

#anim = @animate for i=2:2:360
anim = @animate for i=10:10:360

  if rem(i,20) == 0
    ORCA.restart_server()
  end
  plot!(camera=(i,30))
end
gif(anim, "leastsq-quadratic-samples-log.gif")

using LinearAlgebra
plotlyjs()
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)/size(A,1)
ezsurf!(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot!(x, y, Z, st=:surface, alpha=0.4)
end

## Show the same thing with a log10
quadratic_leastsq(A,b) = x-> log10(0.5*norm(A*x - b)^2)/size(A,1)

Random.seed!(0)
nsamp = 10
k = 25

plot(framestyle=:none,colorbar=false,size=(300,300),dpi=300)
plot!(colorbar=false,size=(900,900))

Random.seed!(0)
for j=1:k
  subset = randperm(size(A,1))[1:nsamp] # choose 25 of 200 samples.
  ezsurf!(-3:0.05:3,-50:0.25:250,
    quadratic_leastsq(A[subset,:],b[subset,:]))
end

anim = @animate for i=2:2:360
#anim = @animate for i=10:10:360
  if rem(i,20) == 0
    ORCA.restart_server()
  end
  plot!(camera=(i,30))
end
gif(anim, "leastsq-quadratic-samples-log.gif")


text:
struct SlicedArray{T<:AbstractArray,N,M,P<:AbstractArray,A<:TupleN{TypedBool},C} <: AbstractArray{T,N}
    parent::P
    alongs::A
    @inline function SlicedArray{T,N,M,P,A,C}(
        parent::P,
        alongs::A,
    ) where {T<:AbstractArray,N,M,P<:AbstractArray,A<:TupleN{TypedBool},C}
        # TODO check type parameters
        new(parent, alongs)
    end
end

@inline function SlicedArray(
    parent::P,
    alongs::A,
    paxes::Anys{L},
    outaxes::Anys{N},
    inaxes::Anys{M},
) where {L,P<:AbstractArray{<:Any,L},A<:NTuple{L,TypedBool},N,M}
    I = tuple_setindex(sliceaxes(parent), tuple_map(first, outaxes), invert(alongs))
    T = viewtype(parent, I)
    SlicedArray{T,N,M,P,A,iscontiguous(alongs)}(parent, alongs)
end

@inline function SlicedArray(parent::AbstractArray{<:Any,L}, alongs::NTuple{L,TypedBool}) where {L}
    paxes = sliceaxes(parent)
    inaxes = tuple_getindex(paxes, alongs)
    outaxes = tuple_getindex(paxes, tuple_map(!, alongs))
    SlicedArray(parent, alongs, paxes, outaxes, inaxes)
end

####
#### Core Array Interface
####

@inline Base.axes(S::SlicedArray) = tuple_getindex(axes(S.parent), invert(S.alongs))

@inline Base.size(S::SlicedArray) = tuple_getindex(size(S.parent), invert(S.alongs))

@propagate_inbounds function Base.getindex(S::SlicedArray{<:Any,N}, I::Vararg{Int,N}) where {N}
    view(S.parent, mergeindices(S, I)...)
end

@propagate_inbounds function Base.setindex!(S::SlicedArray{<:Any,N}, v, I::Vararg{Int,N}) where {N}
    S.parent[mergeindices(S, I)...] = v
    return S
end

# eltype(S) <: AbstractArray{T,0}, but eltype(S.parent) is just T, so we need to unwrap
# zero-dimensional elements.
@propagate_inbounds function Base.setindex!(
    S::SlicedArray{<:Any,N,0},
    v::AbstractArray{<:Any,0},
    I::Vararg{Int,N},
) where {N}
    S.parent[mergeindices(S, I)...] = v[]
    return S
end

Base.IndexStyle(::Type{<:SlicedArray}) = IndexCartesian()
Base.IndexStyle(::Type{<:SlicedArray{<:Any,1}}) = IndexLinear()

@inline function mergeindices(S::SlicedArray{<:Any,N}, I::Anys{N}) where {N}
    tuple_setindex(sliceaxes(S.parent), I, invert(S.alongs))
end

function Base.resize!(S::SlicedArray{<:Any,N}, dims::NTuple{N,Integer}) where {N}
    resize!(S.parent, tuple_setindex(size(S.parent), dims, invert(S.alongs)))
    return S
end
Base.resize!(S::SlicedArray{<:Any,N}, dims::Vararg{Integer,N}) where {N} = resize!(S, dims)

function Base.sizehint!(S::SlicedArray{<:Any,N}, dims::NTuple{N,Integer}) where {N}
    sizehint!(S.parent, tuple_setindex(size(S.parent), dims, invert(S.alongs)))
    return S
end
Base.sizehint!(S::SlicedArray{<:Any,N}, dims::Vararg{Integer,N}) where {N} = sizehint!(S, dims)

@inline function setindex_shape_check(dest::SlicedArray, szsrc::Dims)
    setindex_shape_check(inner_size(dest), szsrc)
end

####
#### ContiguousSlicedArray
####

# If only the leading dimensions of S.parent are sliced (i.e. S.alongs consists of any number of
# True's followed by any number of False's) then we can optimize some functions, and if N==1
# then functions like push! are well-defined, assuming the underlying parent supports those
# operations (e.g. ElasticArray).

const ContiguousSlicedArray{T,N,M,P,A} = SlicedArray{T,N,M,P,A,true}
const ContiguousSlicedVector{T,M,P,A} = ContiguousSlicedArray{T,1,M,P,A}

function Base.copyto!(
    dest::ContiguousSlicedArray{<:Any,N,M},
    doffs::Integer,
    src::ContiguousSlicedArray{<:Any,N,M},
    soffs::Integer,
    n::Integer,
) where {N,M}
    insize = inner_size(src)
    setindex_shape_check(dest, insize)
    stride = prod(insize)
    di1 = firstindex(dest.parent)
    si1 = firstindex(src.parent)
    doffs_parent = stride * (doffs - di1) + di1
    soffs_parent = stride * (soffs - si1) + si1
    copyto!(dest.parent, doffs_parent, src.parent, soffs_parent, stride * n)
    return dest
end

function Base.copyto!(
    dest::ContiguousSlicedArray{<:Any,N,M},
    src::ContiguousSlicedArray{<:Any,N,M},
) where {N,M}
    setindex_shape_check(dest, inner_size(src))
    copyto!(dest.parent, src.parent)
    return dest
end

function Base.append!(S::ContiguousSlicedVector, iter::ContiguousSlicedVector)
    setindex_shape_check(S, inner_size(iter))
    append!(S.parent, iter.parent)
    return S
end

# TODO can probably speed this up by doing resize + setindex like in Base
function Base.append!(S::ContiguousSlicedVector, iter)
    for el in iter
        push!(S, el)
    end
    return S
end

function Base.prepend!(S::ContiguousSlicedVector, iter::ContiguousSlicedVector)
    setindex_shape_check(S, inner_size(iter))
    prepend!(S.parent, iter.parent)
    return S
end

# TODO can probably speed this up by doing resize + setindex like in Base
function Base.prepend!(S::ContiguousSlicedVector, iter)
    for el in iter
        pushfirst!(S, el)
    end
    return S
end

function Base.pop!(S::ContiguousSlicedVector)
    isempty(S) && argerror("array must be non-empty")
    # We need to copy the last element as the view could be invalidated by resizing S.parent,
    # but Base.copy(::SubArray{T,0}) unpacks the 0-dimensional view and an element of type T,
    # so we wrap it with a 0-dimensional array
    item = _maybe_wrap(eltype(S.parent), copy(last(S)))
    resize!(S, length(S) - 1)
    return item
end

# TODO need something like resizebeg! for ElasticArray
#function Base.popfirst!(S::ContiguousSlicedVector)
#    isempty(S) && argerror("array must be non-empty"))
#    item = _maybe_wrap(eltype(S.parent), copy(first(S)))
#    resizebeg!(S.parent, length(S) - 1)
#    return item
#end

@inline _maybe_wrap(::Type{T}, x::AbstractArray{T}) where {T} = x
@inline _maybe_wrap(::Type{T}, x::T) where {T} = fill(x)

function Base.push!(S::ContiguousSlicedVector{<:Any,M}, item::AbstractArray{<:Any,M}) where {M}
    setindex_shape_check(S, size(item))
    append!(S.parent, item)
    return S
end

function Base.pushfirst!(S::ContiguousSlicedVector{<:Any,M}, item::AbstractArray{<:Any,M}) where {M}
    setindex_shape_check(S, size(item))
    prepend!(S.parent, item)
    return S
end

Base.empty!(S::ContiguousSlicedVector) = resize!(S, 0)

flatten(S::ContiguousSlicedArray) = S.parent

#####
##### Misc
#####

Base.parent(S::SlicedArray) = S.parent

Base.dataids(S::SlicedArray) = Base.dataids(S.parent)

Base.copy(S::SlicedArray) = SlicedArray(copy(S.parent), S.alongs)

function Base.showarg(io::IO, S::SlicedArray, toplevel)
    print(io, "slice(")
    Base.showarg(io, S.parent, false) # TODO try ture
    if length(S.alongs) > 0
        print(io, ", ", join(map(a -> a === True() ? ':' : '*', S.alongs), ", "))
    end
    print(io, ')')
    if toplevel
        print(io, " with ", Base.dims2string(inner_size(S)), " elements of ", inner_eltype(S))
    end
    return nothing
end

function mapslices(f, A::AbstractArray; dims)
    Asliced = slice(A, dims...)
    S = _alloc_mapslices(Asliced, f(_unsafe_getindex_unwrapped(Asliced, firstindex(Asliced))))
    for I in eachindex(Asliced, S) # TODO start from second
        _unsafe_setindex!(S, f(_unsafe_getindex_unwrapped(Asliced, I)), I)
    end
    return S.parent
end

@inline function _alloc_mapslices(S::SlicedArray{<:Any,N,M}, b1) where {N,M}
    insize = tuple_map(unsafe_length, _reshape_axes(axes(b1), Val(M)))
    psize = tuple_setindex(size(S.parent), insize, S.alongs)
    SlicedArray(Array{eltype(b1),M + N}(undef, psize...), S.alongs)
end
@inline _reshape_axes(axes::Tuple, ::Val{N}) where {N} = Base.rdims(Val(N), axes)
@inline _reshape_axes(axes::NTuple{N,Any}, ::Val{N}) where {N} = axes

@inline _unsafe_getindex_unwrapped(S::SlicedArray{<:Any,<:Any,0}, I) = @inbounds S[I][]
@inline _unsafe_getindex_unwrapped(S::SlicedArray, I) = @inbounds S[I]

@inline _unsafe_setindex!(S::SlicedArray, v, i::Int) = @inbounds S[i] .= v
@inline _unsafe_setindex!(S::SlicedArray, v, I::CartesianIndex) = @inbounds S[Tuple(I)...] .= v
@inline _unsafe_setindex!(S::SlicedArray, v::AbstractArray, I::CartesianIndex) = @inbounds S[I] = v
@inline _unsafe_setindex!(S::SlicedArray, v::AbstractArray, i::Int) = @inbounds S[i] = v

#####
##### Extra
#####

@inline inner_axes(S::SlicedArray) = tuple_getindex(axes(S.parent), S.alongs)

@inline inner_size(S::SlicedArray) = tuple_getindex(size(S.parent), S.alongs)

"""
    slice(A::AbstractArray, alongs...)
    slice(A::AbstractArray, alongs)

Return an array whose elements are views into `A` along the dimensions `alongs`.
`alongs` can be specified in the following ways:
- `:` or `*` representing sliced or indexed dimensions, respectively.
- Integers corresponding to the sliced dimensions.

See also: [`align`](@ref), [`flatten`](@ref).

# Examples

```jldoctest
julia> A = reshape(collect(1:8), (2,2,2))
2×2×2 Array{Int64,3}:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> S = slice(A, :, *, :)
2-element slice(::Array{Int64,3}, :, *, :) with 2×2 elements of Int64:
 [1 5; 2 6]
 [3 7; 4 8]

julia> S[1] == view(A, :, 1, :)
true

julia> slice(A, :, *, :) == slice(A, 1, 3)
true
```
"""
@inline function slice(A::AbstractArray{<:Any,L}, alongs...) where {L}
    _slice(A, to_alongs(Val(L), alongs...))
end
@inline function slice(A::AbstractArray{<:Any,L}, alongs::Tuple) where {L}
    _slice(A, to_alongs(Val(L), alongs))
end

@inline function _slice(A::AbstractArray{<:Any,L}, alongs::NTuple{L,TypedBool}) where {L}
    SlicedArray(A, alongs)
end
# When taking a slice of a slice we can avoid the extra indirection by merging alongs with S.alongs
@inline function _slice(S::SlicedArray{<:Any,N}, alongs::NTuple{N,TypedBool}) where {N}
    palongs = tuple_setindex(S.alongs, alongs, tuple_map(!, S.alongs))
    SlicedArray(S.parent, palongs)
end

"""
    align(A::AbstractArray{<:AbstractArray{V,M},N}, alongs...)
    align(A::AbstractArray{<:AbstractArray{V,M},N}, alongs)

The inverse of `slice`. Returns an `M+N`-dimension array where `alongs` specifies the dimensions
taken up by the inner arrays. If `A` is not a subtype of `AbstractArray{<:AbstractArray{V,M},N}`,
the return value is `A` itself.

See also: [`flatten`](@ref), [`slice`](@ref).

# Examples

```jldoctest
julia> A = [[1 2; 3 4], [5 6; 7 8]]
2-element Array{Array{Int64,2},1}:
 [1 2; 3 4]
 [5 6; 7 8]

julia> align(A, :, :, *)
2×2×2 PermutedDimsArray(flatten(::Array{Array{Int64,2},1}), (1, 2, 3)) with eltype Int64:
[:, :, 1] =
 1  2
 3  4

[:, :, 2] =
 5  6
 7  8

julia> align(A, :, *, :)
2×2×2 PermutedDimsArray(flatten(::Array{Array{Int64,2},1}), (1, 3, 2)) with eltype Int64:
[:, :, 1] =
 1  5
 3  7

[:, :, 2] =
 2  6
 4  8

julia> align(A, *, :, :) == align(A, 2, 3)
true

julia> align(slice(A, 2, 3), 2, 3) == A
true
```
"""
@inline function align(A::NestedArray{<:Any,M,N}, alongs...) where {M,N}
    _align(A, to_alongs(Val(M+N), alongs...))
end
@inline function align(A::NestedArray{<:Any,M,N}, alongs::Tuple) where {M,N}
    _align(A, to_alongs(Val(M+N), alongs))
end

# align(A, alongs) is the inverse of slice(A, alongs) i.e. align(slice(A, alongs), alongs) === A
function _align(S::SlicedArray{<:Any,<:Any,<:Any,<:Any,A}, alongs::A) where {A<:TupleN{TypedBool}}
    S.parent
end
@inline function _align(A::NestedArray{<:Any,M,N}, alongs::NTuple{L,TypedBool}) where {M,N,L}
    (M+N == L && M == length(tuple_getindex(alongs, alongs))) || throw_invalid_alongs(M, N, alongs)
    return iscontiguous(alongs) ? flatten(A) : PermutedDimsArray(flatten(A), _permdims(alongs))
end

@pure function _permdims(alongs::NTuple{L,TypedBool}) where {L}
    dims = ntuple(identity, Val(L))
    return (tuple_getindex(dims, alongs)..., tuple_getindex(dims, invert(alongs))...)
end

#####
##### 3rd Party
#####

@inline function UnsafeArrays.unsafe_uview(S::SlicedArray)
    SlicedArray(UnsafeArrays.unsafe_uview(S.parent), S.alongs)
end

function Adapt.adapt_structure(T, S::SlicedArray)
    SlicedArray(adapt(T, S.parent), S.alongs)
end


text:
using OrbitalInsertion, InterpHalo, BenchmarkTools

# Create the Halo orbit with InterpHalo
ORB = map(x -> intH(x,0.5,2), range(0.0,stop=1.0,length=201) )
orbita = zeros(length(ORB),6)
for i in 1:size(orbita,1)
    orbita[i,:] = ORB[i]
end
# Select one point of the orbit
y0 = point(1,orbita)[1]

# Create the problem relevant magnitudes
sys         = createSystem(EMsystem()) # when EMsystem is used it gives back the properties of the Earth-Moon system
prob_info   = createConstants(sys)  # which are now used to create the problem information object needed for the optimization

firstGuess = [0.05,-0.05,0.1]
par = createPar(prob_info;  tspan = -2.8, hgoal = 300, hmmin = 100 ,
                            sType = "global", points = 30, x0 = firstGuess,
                            maxt0 = -1.0, rest = "basic", lim = 3000,
                            cTol = 1e-5, tRest = "mixed", alg = "lowOrder",
                            atol = 1e-9, rtol = 1e-9, retorno = false, detMin = false, y = y0, optTol = 1e-8, neval = 450)

gsol = search(par,prob_info)
# @time search(par,prob_info) # test the time to compute the global search

fsol, sFsol = filterGlobalSol(gsol, par,prob_info)

plot = plotSols(sFsol, par, prob_info, orbita,true)

# Make local search with the best of the global solutions
X0 = sFsol[2][2]

par = createPar(prob_info;  tspan = -2.8, hgoal = 300, hmmin = 100 ,
                            sType = "local", x0 = X0,
                            maxt0 = -1.0, rest = "basic", lim = 3000,
                            cTol = 1e-10, tRest = "mixed", alg = "lowOrder",
                            atol = 1e-14, rtol = 1e-14, retorno = false, detMin = false, y = y0, optTol = 1e-12, neval = 1500)

finalSolution = search(par,prob_info)
plot = plotSols(finalSolution, par, prob_info, orbita,true)