text:
@testset "Renewable data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    info_message = "The data doesn't include devices of type HydroEnergyReservoir, consider changing the device models"
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    @test_logs (:info, info_message) match_mode = :any mock_construct_device!(
        model,
        device_model,
    )
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys14)
    @test_logs (:info, info_message) match_mode = :any mock_construct_device!(
        model,
        device_model,
    )
end

###################################
###### FIXED OUTPUT TESTS #########
###################################

@testset "Hydro DCPLossLess FixedOutput" begin
    device_model = DeviceModel(HydroDispatch, FixedOutput)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro DCPLossLess HydroEnergyReservoir with FixedOutput formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, FixedOutput)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

###################################
### RUN OF RIVER DISPATCH TESTS ###
###################################

@testset "Hydro DCPLossLess HydroDispatch with HydroDispatchRunOfRiver formulations" begin
    device_model = DeviceModel(HydroDispatch, HydroDispatchRunOfRiver)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 24, 0, 48, 24, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPPowerModel HydroDispatch with HydroDispatchRunOfRiver formulations" begin
    device_model = DeviceModel(HydroDispatch, HydroDispatchRunOfRiver)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 48, 0, 72, 48, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro DCPLossLess HydroEnergyReservoir with HydroDispatchRunOfRiver formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 24, 0, 48, 24, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPPowerModel HydroEnergyReservoir with HydroDispatchRunOfRiver formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 48, 0, 72, 48, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

###################################
#### RUN OF RIVER COMMIT TESTS ####
###################################

@testset "Hydro DCPLossLess HydroDispatch with HydroCommitmentRunOfRiver formulations" begin
    device_model = DeviceModel(HydroDispatch, HydroCommitmentRunOfRiver)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 48, 0, 48, 24, 0, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPPowerModel HydroDispatch with HydroCommitmentRunOfRiver formulations" begin
    device_model = DeviceModel(HydroDispatch, HydroCommitmentRunOfRiver)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 72, 0, 72, 48, 0, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro DCPLossLess HydroEnergyReservoir with HydroCommitmentRunOfRiver formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentRunOfRiver)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 48, 0, 48, 24, 0, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPPowerModel HydroEnergyReservoir with HydroCommitmentRunOfRiver formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentRunOfRiver)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 72, 0, 72, 48, 0, true)
    psi_checkobjfun_test(model, GAEVF)
end

#########################################
#### RESERVOIR BUDGET DISPATCH TESTS ####
#########################################

@testset "Hydro DCPLossLess HydroEnergyReservoir with HydroDispatchReservoirBudget Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirBudget)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 24, 0, 25, 24, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPPowerModel HydroEnergyReservoir with HydroDispatchReservoirBudget Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirBudget)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 48, 0, 49, 48, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

#########################################
#### PUMPED STORAGE DISPATCH TESTS ####
#########################################

@testset "Hydro DCPLossLess HydroPumpedStorage with HydroDispatchPumpedStorage Formulations" begin
    device_model = DeviceModel(
        HydroPumpedStorage,
        HydroDispatchPumpedStorage;
        attributes=Dict{String, Any}("reservation" => false),
    )
    c_sys5_phes_ed = PSB.build_system(PSITestSystems, "c_sys5_phes_ed")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_phes_ed)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 60, 0, 24, 24, 24, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro DCPLossLess HydroPumpedStorage with HydroDispatchPumpedStorage with Reservation Formulations" begin
    device_model = DeviceModel(HydroPumpedStorage, HydroDispatchPumpedStorage)
    c_sys5_phes_ed = PSB.build_system(PSITestSystems, "c_sys5_phes_ed")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_phes_ed)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 72, 0, 24, 24, 24, true)
    psi_checkobjfun_test(model, GAEVF)
end

#########################################
### RESERVOIR BUDGET COMMITMENT TESTS ###
#########################################

@testset "Hydro DCPLossLess HydroEnergyReservoir with HydroCommitmentReservoirBudget Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentReservoirBudget)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 48, 0, 25, 24, 0, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPPowerModel HydroEnergyReservoir with HydroCommitmentReservoirBudget Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentReservoirBudget)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 72, 0, 49, 48, 0, true)
    psi_checkobjfun_test(model, GAEVF)
end

#########################################
### RESERVOIR STORAGE DISPATCH TESTS ####
#########################################

@testset "Hydro DCPLossLess HydroEnergyReservoir with HydroDispatchReservoirStorage Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 24, 24, 48, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPLossLess HydroEnergyReservoir with HydroDispatchReservoirStorage Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 144, 0, 48, 48, 48, false)
    psi_checkobjfun_test(model, GAEVF)
end

#########################################
### RESERVOIR STORAGE COMMITMENT TESTS ##
#########################################

@testset "Hydro DCPLossLess HydroEnergyReservoir with HydroCommitmentReservoirStorage Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentReservoirStorage)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 144, 0, 24, 24, 48, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hydro ACPLossLess HydroEnergyReservoir with HydroCommitmentReservoirStorage Formulations" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentReservoirStorage)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_hyd)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 168, 0, 48, 48, 48, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Solving ED Hydro System using Dispatch Run of River" begin
    sys = PSB.build_system(PSITestSystems, "c_sys5_hy")
    networks = [ACPPowerModel, DCPPowerModel]

    test_results = Dict{Any, Float64}(ACPPowerModel => 177526.0, DCPPowerModel => 175521.0)

    for net in networks
        @testset "HydroRoR ED model $(net)" begin
            template = get_thermal_dispatch_template_network(net)
            set_device_model!(template, HydroDispatch, HydroDispatchRunOfRiver)
            ED = DecisionModel(
                EconomicDispatchProblem,
                template,
                sys;
                optimizer=ipopt_optimizer,
            )
            @test build!(ED; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
            psi_checksolve_test(
                ED,
                [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
                test_results[net],
                1000,
            )
        end
    end
end

@testset "Solving ED Hydro System using Commitment Run of River" begin
    sys = PSB.build_system(PSITestSystems, "c_sys5_hy")
    net = DCPPowerModel

    template = get_thermal_dispatch_template_network(net)
    set_device_model!(template, HydroDispatch, HydroCommitmentRunOfRiver)

    @testset "HydroRoR ED model $(net)" begin
        ED = DecisionModel(UnitCommitmentProblem, template, sys; optimizer=GLPK_optimizer)
        @test build!(ED; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_checksolve_test(ED, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], 175521.0, 1000)
    end
end

@testset "Solving ED Hydro System using Dispatch with Reservoir" begin
    systems = [
        PSB.build_system(PSITestSystems, "c_sys5_hyd"),
        PSB.build_system(PSITestSystems, "c_sys5_hyd_ems"),
    ]
    networks = [ACPPowerModel, DCPPowerModel]
    models = [HydroDispatchReservoirBudget, HydroDispatchReservoirStorage]
    test_results = Dict{Any, Float64}(
        (ACPPowerModel, HydroDispatchReservoirBudget) => 33423.0,
        (DCPPowerModel, HydroDispatchReservoirBudget) => 33042.0,
        (ACPPowerModel, HydroDispatchReservoirStorage) => 232497.0,
        (DCPPowerModel, HydroDispatchReservoirStorage) => 230153.0,
    )

    for net in networks, (mod, sys) in zip(models, systems)
        @testset "$(mod) ED model on $(net)" begin
            template = get_thermal_dispatch_template_network(net)
            set_device_model!(template, HydroEnergyReservoir, mod)

            ED = DecisionModel(
                EconomicDispatchProblem,
                template,
                sys;
                optimizer=ipopt_optimizer,
            )
            @test build!(ED; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
            psi_checksolve_test(
                ED,
                [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
                test_results[(net, mod)],
                10000,
            )
        end
    end
end

@testset "Solving ED Hydro System using Commitment with Reservoir" begin
    systems = [
        PSB.build_system(PSITestSystems, "c_sys5_hyd"),
        PSB.build_system(PSITestSystems, "c_sys5_hyd_ems"),
    ]
    net = DCPPowerModel
    models = [HydroCommitmentReservoirBudget, HydroCommitmentReservoirStorage]
    test_results = Dict{Any, Float64}(
        HydroCommitmentReservoirBudget => 33042.0,
        HydroCommitmentReservoirStorage => 230153.0,
    )

    for (mod, sys) in zip(models, systems)
        @testset "$(mod) ED model on $(net) and use_parameters = true" begin
            template = get_thermal_dispatch_template_network(net)
            set_device_model!(template, HydroEnergyReservoir, mod)

            ED = DecisionModel(
                UnitCommitmentProblem,
                template,
                sys;
                optimizer=GLPK_optimizer,
            )
            @test build!(ED; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
            psi_checksolve_test(
                ED,
                [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
                test_results[mod],
                10000,
            )
        end
    end
end

@testset "HydroEnergyReservoir with HydroDispatchReservoirStorage Formulations (energy target - cases 1b-2b)" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "hydro_test_case_b_sys")

    model = DecisionModel(
        EconomicDispatchProblem,
        template,
        c_sys5_hyd;
        optimizer=HiGHS_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 15, 0, 3, 3, 9, false)
    psi_checksolve_test(model, [MOI.OPTIMAL], 5621.0, 10.0)
end

@testset "HydroEnergyReservoir with HydroDispatchReservoirStorage Formulations (energy target - cases 1c-2c)" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "hydro_test_case_c_sys")

    model = DecisionModel(
        EconomicDispatchProblem,
        template,
        c_sys5_hyd;
        optimizer=HiGHS_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 15, 0, 3, 3, 9, false)
    psi_checksolve_test(model, [MOI.OPTIMAL], 21.0)
end

@testset "HydroEnergyReservoir with HydroDispatchReservoirStorage Formulations (energy target - cases 1d-2d)" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "hydro_test_case_d_sys")

    model = DecisionModel(
        EconomicDispatchProblem,
        template,
        c_sys5_hyd;
        optimizer=HiGHS_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 15, 0, 3, 3, 9, false)
    psi_checksolve_test(model, [MOI.OPTIMAL], -5429.0, 10.0)
end

@testset "HydroEnergyReservoir with HydroDispatchReservoirStorage Formulations (energy target - cases 1e-2e)" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "hydro_test_case_e_sys")

    model = DecisionModel(
        EconomicDispatchProblem,
        template,
        c_sys5_hyd;
        optimizer=HiGHS_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 15, 0, 3, 3, 9, false)
    psi_checksolve_test(model, [MOI.OPTIMAL], 21.0, 10.0)
end

@testset "HydroEnergyReservoir with HydroDispatchReservoirStorage Formulations (energy target - cases 1f-2f)" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5_hyd = PSB.build_system(PSITestSystems, "hydro_test_case_f_sys")

    model = DecisionModel(
        EconomicDispatchProblem,
        template,
        c_sys5_hyd;
        optimizer=HiGHS_optimizer,
    )
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 15, 0, 3, 3, 9, false)
    psi_checksolve_test(model, [MOI.OPTIMAL], -17179.0)
end

### Feedforward Test ###

@testset "Test SemiContinuousFeedforward to HydroDispatch with HydroCommitmentRunOfRiver model" begin
    device_model = DeviceModel(HydroDispatch, HydroCommitmentRunOfRiver)
    ff_sc = SemiContinuousFeedforward(
        component_type=HydroDispatch,
        source=OnVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_sc)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 48, 0, 72, 48, 0, true)
end

@testset "Test UpperBoundFeedforward to HydroDispatch with HydroCommitmentRunOfRiver model" begin
    device_model = DeviceModel(HydroDispatch, HydroCommitmentRunOfRiver)
    ff_ub = UpperBoundFeedforward(
        component_type=HydroDispatch,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 48, 0, 72, 24, 0, true)
end

@testset "Test LowerBoundFeedforward to HydroDispatch with HydroCommitmentRunOfRiver model" begin
    device_model = DeviceModel(HydroDispatch, HydroCommitmentRunOfRiver)
    ff_ub = LowerBoundFeedforward(
        component_type=HydroDispatch,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 48, 0, 48, 48, 0, true)
end

@testset "Test UpperBoundFeedforward to HydroDispatch with HydroDispatchRunOfRiver model" begin
    device_model = DeviceModel(HydroDispatch, HydroDispatchRunOfRiver)

    ff_ub = UpperBoundFeedforward(
        component_type=HydroDispatch,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 24, 0, 72, 24, 0, false)
end

@testset "Test LowerBoundFeedforward to HydroDispatch with HydroDispatchRunOfRiver model" begin
    device_model = DeviceModel(HydroDispatch, HydroDispatchRunOfRiver)

    ff_ub = LowerBoundFeedforward(
        component_type=HydroDispatch,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hy")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 24, 0, 48, 48, 0, false)
end

@testset "Test EnergyLimitFeedforward to HydroEnergyReservoir with HydroDispatchReservoirBudget model" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirBudget)
    ff_il = EnergyLimitFeedforward(
        component_type=HydroEnergyReservoir,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
        number_of_periods=12,
    )

    PSI.attach_feedforward!(device_model, ff_il)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hyd")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 24, 0, 27, 24, 0, false)
end

@testset "Test EnergyLimitFeedforward to HydroEnergyReservoir with HydroDispatchReservoirStorage model" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage)
    ff_il = EnergyLimitFeedforward(
        component_type=HydroEnergyReservoir,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
        number_of_periods=12,
    )

    PSI.attach_feedforward!(device_model, ff_il)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 26, 24, 48, false)
end

@testset "Test LowerBoundFeedforward to HydroEnergyReservoir with HydroDispatchReservoirStorage model" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage)
    ff_il = LowerBoundFeedforward(
        component_type=HydroEnergyReservoir,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_il)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 24, 48, 48, false)
end

@testset "Test EnergyLimitFeedforward to HydroEnergyReservoir with HydroCommitmentReservoirStorage model" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentReservoirStorage)
    ff_il = EnergyLimitFeedforward(
        component_type=HydroEnergyReservoir,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
        number_of_periods=12,
    )
    PSI.attach_feedforward!(device_model, ff_il)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 144, 0, 26, 24, 48, true)
end

@testset "Test SemiContinuousFeedforward to HydroEnergyReservoir with HydroCommitmentReservoirStorage model" begin
    device_model = DeviceModel(HydroEnergyReservoir, HydroCommitmentReservoirStorage)
    ff_sc = SemiContinuousFeedforward(
        component_type=HydroEnergyReservoir,
        source=OnVariable,
        affected_values=[ActivePowerVariable],
    )
    PSI.attach_feedforward!(device_model, ff_sc)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_hyd_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 144, 0, 48, 48, 48, true)
end

@testset "Test EnergyLimitFeedforward to HydroEnergyReservoir models" begin
    device_model = DeviceModel(HydroPumpedStorage, HydroDispatchPumpedStorage)

    ff_il = EnergyLimitFeedforward(
        component_type=HydroPumpedStorage,
        source=ActivePowerOutVariable,
        affected_values=[ActivePowerOutVariable],
        number_of_periods=12,
    )

    PSI.attach_feedforward!(device_model, ff_il)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_phes_ed")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 72, 0, 25, 24, 24, true)
end

@testset "Test EnergyTargetFeedforward to HydroEnergyReservoir models" begin
    device_model = DeviceModel(HydroPumpedStorage, HydroDispatchPumpedStorage)

    ff_up = EnergyTargetFeedforward(
        component_type=HydroPumpedStorage,
        source=EnergyVariableUp,
        affected_values=[EnergyVariableUp],
        target_period=12,
        penalty_cost=1e4,
    )

    PSI.attach_feedforward!(device_model, ff_up)
    c_sys5_hy = PSB.build_system(PSITestSystems, "c_sys5_phes_ed")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_hy)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 84, 0, 24, 25, 24, true)
end


text:
module NoiseReduction

using SFRDSP
using SampledSignals
using DSP
using TensorCast
using Statistics
using TSVD
using LinearAlgebra

export spec_sub, wiener, lowrank

function _spec_sub(X, fr, smoothing, α, β)
    # implements spectral subtraction from (Martin 1994)
    X2 = abs2.(X)
    npow = noise_psd(X, smoothing, fr)
    # α is an "oversubtract" parameter to compensate for the underestimated
    # noise power and is typically in the range of 2-6
    # β is a "noise floor" parameter that will limit how low this will push result,
    # and is usually in the range  0 < β << 1
    # taking the minimum of X2 and the noise floor here is a slight tweak on the
    # original in (berouti 1979), because if the signal happens to be below the noise
    # noise floor the given equation would end up with a bin with greater energy than
    # the mixture. This is particularly problematic when you're using this to estimate
    # the wiener filter.
    X2_sub = @. max(min(X2, β*npow), X2 - α * npow)
    @. sqrt(X2_sub) * exp(im*angle(X))
end

"""
Perform spectral subtraction by first estimating the noise power spectral density, then
subtracting the noise energy from the mixture energy.

`audio`            - a SampleBuf
`nfft`             - The FFT size of the STFT
`hop`              - The hop size of the STFT
`noise_smoothing`  - The window size (in seconds) used to smooth each band before taking the
                     minimum
`oversub`          - Apply this multiplier to the noise estimate before subtracting. Typically
                     in the range of 1-6.
`specfloor`        - Use this (multiplied by the noise power) as a lower bound on each band
                     after subtraction. This is used to limit "musical noise". Typically in
                     the range of 0.005 to 0.1.
"""
function spec_sub(audio;
                  nfft=1024, hop=nfft÷2,
                  noise_smoothing=0.25, oversub=4, specfloor=0.01)
    fs = samplerate(audio)
    X = stft2(audio, nfft, hop; window=cosine)
    X_ss = _spec_sub(X, fs/hop, noise_smoothing, oversub, specfloor)
    # N_ss = X - X_ss
    x_ss = SampleBuf(istft2(X_ss, nfft, hop; window=cosine)[1:length(audio)], fs)
    # n_ss = SampleBuf(istft2(N_ss, nfft, hop; window=cosine)[1:length(audio)], fs)
    n_ss = audio-x_ss
    x_ss, n_ss
end

"""
Perform noise reduction by Wiener Filtering

`audio`            - a SampleBuf
`nfft`             - The FFT size of the STFT
`hop`              - The hop size of the STFT
`noise_smoothing`  - The window size (in seconds) used to smooth each band before taking the
                     minimum
`oversub`          - Apply this multiplier to the noise estimate before subtracting. Typically
                     in the range of 1-6.
`specfloor`        - Use this (multiplied by the noise power) as a lower bound on each band
                     after subtraction. This is used to limit "musical noise". Typically in
                     the range of 0.005 to 0.1.
`wf_smoothing`     - 2D Smoothing applied to the wiener filter before applying (in
                     samples). Set to `1` to disable smoothing.
"""

function wiener(audio; nfft=1024, hop=nfft÷2,
                noise_smoothing=0.25,
                oversub=4, specfloor=0.01, wf_smoothing=11)
    @assert isodd(wf_smoothing)
    fs = samplerate(audio)
    wf_L = wf_smoothing
    X = stft2(audio, nfft, hop; window=cosine)
    # note that we could do this more efficiently in a single step rather than
    # doing the spectral subtraction and then computing the wiener filter
    # coefficient. Not sure if we can incorporate the α and β params though
    X_ss = _spec_sub(X, fs/hop, noise_smoothing, oversub, specfloor)

    # compute and smooth the wiener filter
    wf = @. abs2(X_ss) / abs2(X)
    win = gaussian((wf_L, wf_L), (0.15, 0.15))
    win ./= sum(win)
    wf_sm = conv(win, wf)[wf_L÷2+1:end-wf_L÷2, wf_L÷2+1:end-wf_L÷2]

    # apply the wiener filter
    X_wf = X .* wf_sm

    x_wf = SampleBuf(istft2(X_wf, nfft, hop; window=cosine)[1:length(audio)], fs)
    x_wf, audio-x_wf
end

"""
Performs a low-rank projection across time and channels in each frequency band. Assumes
the input audio has been pre-aligned.

Returns a `Vector` of de-noised signals, as well as a `Vector` of the noise channels,
which are each channel projected into the noise subspace.
"""
function lowrank(audio, rank=1; nfft=1024, hop=nfft÷2,
                 tol=1e-5, niters=20,
                 weight_smoothing=31,
                 wiener=true, wiener_alpha=2, l1_weight=false,
                 normalize=:none,
                 weight_pow=0.5,
                 noise_rank=length(audio)-rank)
    fs = first(unique(samplerate.(audio)))
    @assert(all(samplerate.(audio) .== fs))

    @cast Xs[k, n, c] := stft2(audio[c], nfft, hop; window=cosine, demod=true)[k, n]
    @cast Npsd[k, c] := noise_psd(Xs[:, :, c], 0.25, fs/hop)[k]
    K, N, C = size(Xs)

    Xs_wf = if wiener
        @reduce Xpsd[k, c] := mean(n) abs2(Xs[k, n, c])
        @cast wf[k, c] := max(0.0, Xpsd[k, c] - wiener_alpha*Npsd[k, c]) / Xpsd[k, c]
        @cast [k, n, c] := Xs[k, n, c] * wf[k,c]
    else
        Xs
    end

    # normalize each channel so the louder ones don't dominate the SVD.
    # also flip the dimensions here for better performance downstream
    Xs_norm = if normalize == :channel
        @reduce Xpow[c] := mean(k,n) max(1e-8, abs2(Xs_wf[k, n, c]))
        @cast [c, n, k] := Xs_wf[k, n, c] / Xpow[c]
    elseif normalize == :band
        @reduce Xpow[c, k] := mean(n) max(1e-8, abs2(Xs_wf[k, n, c]))
        @cast [c, n, k] := Xs_wf[k, n, c] / Xpow[c, k]
    elseif normalize == :none
        permutedims(Xs_wf, (3,2,1))
    else
        throw(ArgumentError("Unrecognized `normalize` value $normalize"))
    end

    # we smooth the window across frequency so we can accentuate bins that are
    # nearby signal bins above or below - this helps with some of the horizontal
    # banding
    sm_L=weight_smoothing
    win = reshape(gaussian(sm_L, 0.15), 1, :)
    win ./= sum(win)

    # weights for each time-frequency bin
    w = ones(N, K)
    w_last = copy(w)
    dw = Inf

    # some buffers so we don't need to allocate in the loop as much
    Xbuf = similar(Xs_norm, (C, N))
    Xp = similar(Xs_norm, (rank, N))

    # this doesn't necessarily converge monotonically, but hopefully the change
    # gets pretty small
    # if the weight power is 0 then we don't need to iterate, because we're not
    # refining the weights
    while weight_pow > 0 && dw > tol && niters > 0
        for k in 1:K
            # we're computing a low-rank approximation, and updating the weights based
            # on how much of each time point is in our signal subspace.
            @cast Xbuf[c, n] = Xs_norm[c, n, $k] * w[n, $k]
            U, S, V = tsvd(Xbuf, rank)
            @cast Xbuf[c, n] = Xs_norm[c, n, $k]
            mul!(Xp, U', Xbuf)
            # we're weighting by the normalized energy in the signal subspace here.
            # It was previously the square-root, or amplitude ratio. Not clear which
            # is better.
            w[:, k] .= vec(sum(abs2, Xp; dims=1) ./
                           sum(abs2, Xbuf; dims=1)) .^ weight_pow
        end
        # smooth the weights across frequency
        w = conv(win, w)[:, sm_L÷2+1:end-sm_L÷2]

        dw = norm(w .- w_last) / length(w)
        @info "weights changed by $dw"
        w_last .= w
        niters -= 1
    end

    Xs_proj = similar(Xs)
    Ns_proj = similar(Xs)
    # perform the final iteration and project back into the multichannel space
    for k in 1:K
        # X = @view Xs_norm[k, :, :]
        @cast Xbuf[c, n] = Xs_norm[c, n, $k] * w[n, $k]
        # U, S, V = tsvd(Xbuf, rank)
        D = svd(Xbuf)
        Ux = D.U[:, 1:rank]
        Un = D.U[:, end-noise_rank+1:end]

        # previously I had the commented-out projection, but now (Mar 9, 2020) I'm
        # thinking that actually I want to project the signal, rather than just
        # reconstructing from the SVD (which includes the weights.)
        # Xk = D.U[:, 1:rank] * D.S[1:rank] * D.Vt[1:rank, :]
        # nr = rank+C-noise_rank
        # Nk = D.U[:, nr:end] * (D.S[nr:end] .* D.Vt[nr:end, :])

        @cast Xbuf[c, n] = Xs_norm[c, n, $k]
        Xk = Ux * Ux' * Xbuf
        Nk = Un * Un' * Xbuf
        @cast Xs_proj[$k, n, c] = Xk[c, n] * w[n, $k]
        @cast Ns_proj[$k, n, c] = Nk[c, n]
    end

    if normalize == :channel
        @cast Xs_proj[k, n, c] = Xs_proj[k, n, c] * Xpow[c]
        @cast Ns_proj[k, n, c] = Ns_proj[k, n, c] * Xpow[c]
    elseif normalize == :band
        @cast Xs_proj[k, n, c] = Xs_proj[k, n, c] * Xpow[c, k]
        @cast Ns_proj[k, n, c] = Ns_proj[k, n, c] * Xpow[c, k]
    end

    @cast target[c] := SampleBuf(istft2(Xs_proj[:, :, c], nfft, hop;
                                        window=cosine, demod=true),
                                 fs)

    # we also add back the noise removed by the wiener filter
    # TODO: this isn't going to be scaled properly if there was normalization
    @cast Ns_proj[k, n, c] = Ns_proj[k, n, c] + (Xs[k, n, c] - Xs_wf[k, n, c])
    @cast noise[c] := SampleBuf(istft2(Ns_proj[:, :, c], nfft, hop;
                                       window=cosine, demod=true),
                                fs)

    target, noise
end

end # module

ext: jl


text:
include("Newton.jl")
import Base: push!, pop!, length
import Base: show, getindex

"""
DeflationOperator. It is used to describe the following situation. Assume you want to solve `F(x)=0` with a Newton algorithm but you want to avoid the process to return some already known solutions ``roots_i``. You can use `DeflationOperator` to define a function `M(u)` used to find, with Newton iterations, the zeros of the following function
``F(u) / Π_i(dot(u - roots_i, u - roots_i)^{p} + shift) := F(u) / M(u)``. The fields of the struct `DeflationOperator` are as follows:
- power `p`
- `dot` function, this function has to be bilinear and symmetric for the linear solver
- shift
- roots
The deflation operator is is ``M(u) = \\frac{1}{\\prod_{i=1}^{n_{roots}}(shift + norm(u-roots_i)^p)}``
where ``nrm(u) = dot(u,u)``.

Given `defOp::DeflationOperator`, one can access its roots as `defOp[n]` as a shortcut for `defOp.roots[n]`. Also, one can add (resp.remove) a new root by using `push!(defOp, newroot)` (resp. `pop!(defOp)`). Finally `length(defOp)` is a shortcut for `length(defOp.roots)`
"""
struct DeflationOperator{T <: Real, Tf, vectype}
	power::T
	dot::Tf
	shift::T
	roots::Vector{vectype}
end

push!(df::DeflationOperator{T, Tf, vectype}, v::vectype) where {T, Tf, vectype} = push!(df.roots, v)
pop!(df::DeflationOperator) = pop!(df.roots)
getindex(df::DeflationOperator, inds...) = getindex(df.roots, inds...)
length(df::DeflationOperator) = length(df.roots)

function show(io::IO, df::DeflationOperator)
	println(io, "Deflation operator with ", length(df.roots)," roots")
end

function (df::DeflationOperator{T, Tf, vectype})(u::vectype) where {T, Tf, vectype}
	nrm  = u -> df.dot(u, u)
	@assert length(df.roots) > 0 "You need to specify some roots for deflation to work"
	# compute u - df.roots[1]
	tmp = copyto!(similar(u), u);	axpy!(T(-1), df.roots[1], tmp)
	out = T(1) / nrm(tmp)^df.power + df.shift
	for ii = 2:length(df.roots)
		copyto!(tmp, u); axpy!(T(-1), df.roots[ii], tmp)
		out *= T(1) / nrm(tmp)^df.power + df.shift
	end
	return out
end

"""
	pb = DeflatedProblem(F, J, M::DeflationOperator)

This creates a problem encoded a deflated problem ``M(u) \\cdot F(u) = 0`` where `M` is a `DeflationOperator` which encodes the penalization term. `J` is the jacobian of `J`. Can be used to call `newton` and `continuation`.
"""
struct DeflatedProblem{T, Tf, vectype, TF, TJ, def <: DeflationOperator{T, Tf, vectype}}
	F::TF
	J::TJ
	M::def
end

"""
Return the deflated function M(u) * F(u) where M(u) ∈ R
"""
function (df::DeflatedProblem{T, Tf, vectype, TF, TJ, def})(u::vectype) where {T, Tf, vectype, TF, TJ, def}
	out = df.F(u)
	rmul!(out, df.M(u))
	return out
end

###################################################################################################
# Implement the Jacobian operator of the Deflated problem

struct DeflatedLinearSolver <: AbstractLinearSolver end

"""
Implement the linear solver for the deflated problem
"""
function (dfl::DeflatedLinearSolver)(J, rhs)
	# the expression of the Functional is now
	# F(u) / Π_i(dot(u - root_i, u - root_i)^power + shift) := F(u) * M(u)
	# the expression of the differential is
	# dF(u)⋅du * M(u) + F(u) dM(u)⋅du

	# the point at which to compute the Jacobian
	u = J[1]
	# deflated Problem composite type
	defPb = J[2]
	linsolve = J[3]
	Fu = defPb.F(u)
	Mu = defPb.M(u)
	Ju = defPb.J(u)

	# linear solve for the deflated problem. We note that Mu ∈ R
	# hence dM(u)⋅du is a scalar. We now solve the following linear problem
	# M(u) * dF(u)⋅h + F(u) dM(u)⋅h = rhs
	h1, h2, _, (it1, it2) = linsolve(Ju, rhs, Fu)

	# We look for the expression of dM(u)⋅h
	# the solution is then h = Mu * h1 - z h2 where z has to be determined
	delta = 1e-8
	# z1 = (defPb.M(u + delta * h1) - Mu)/delta
	tmp = copyto!(similar(u), u); axpy!(delta, h1, tmp)
	z1 = (defPb.M(tmp) - Mu) / delta

	# z2 = (defPb.M(u + delta * h2) - Mu)/delta
	copyto!(tmp, u); axpy!(delta, h2, tmp)
	z2 = (defPb.M(tmp) - Mu) / delta

	z = z1 / (Mu + z2)

	# return (h1 - z * h2) / Mu, true, (it1, it2)
	copyto!(tmp, h1)
	axpy!(-z, h2, tmp)
	rmul!(tmp, 1 / Mu)
	return tmp, true, (it1, it2)
end
####################################################################################################
"""
	function newton(Fhandle, Jhandle, x0::vectype, options:: NewtonPar{T}, defOp::DeflationOperator{T, Tf, vectype}; kwargs...) where {T, Tf, vectype}

This is the deflated version of the Newton Solver for `F(x) = 0` with Jacobian `J`. It penalises the roots saved in `defOp.roots`. The other arguments are as for `newton`. See [`DeflationOperator`](@ref) for more informations.

Simplified calls are provided, for example when `J` is not passed. It then computed with finite differences.

# Output:
- solution:
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(Fhandle, Jhandle, x0::vectype, options::NewtonPar{T, S, E}, defOp::DeflationOperator{T, Tf, vectype}; kwargs...) where {T, Tf, vectype, S, E}
	# we create the new functional
	deflatedPb = DeflatedProblem(Fhandle, Jhandle, defOp)

	# and its jacobian
	Jacdf = (u0, pb::DeflatedProblem, ls) -> (return (u0, pb, ls))

	# Rename some parameters
	opt_def = @set options.linsolver = DeflatedLinearSolver()
	return newton(u -> deflatedPb(u),
				u -> Jacdf(u, deflatedPb, options.linsolver),
				x0,
				opt_def; kwargs...)
end

# simplified call when no Jacobian is given
function newton(Fhandle, x0::vectype, options::NewtonPar{T, S, E}, defOp::DeflationOperator{T, Tf, vectype};kwargs...) where {T, Tf, vectype, S, E}
	Jhandle = u -> PseudoArcLengthContinuation.finiteDifferences(Fhandle, u)
	return newton(Fhandle,  Jhandle,  x0, options,  defOp;kwargs...)
end


text:
import RvSpectMLBase
using Test

@testset "Spectra1DBasic" begin
    using RvSpectMLBase.TheoreticalInstrument
    n1 = TheoreticalInstrument1D()
    lambdas = TheoreticalInstrument.calc_λs(n1)
    num_pixels = length(lambdas)

    @testset "Create Spectra1DBasic" begin
        @test_nowarn Spectra1DBasic(lambdas, ones(num_pixels), ones(num_pixels), n1)
    end
    spec = Spectra1DBasic(lambdas, ones(num_pixels), ones(num_pixels), n1)
end

@testset "Spectra2DBasic" begin
    num_pixels = 100
    num_orders = 10
    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    n2 = TheoreticalInstrument2D(lambda2d)
    @test_nowarn TheoreticalInstrument.demo_generate_spectrum_line(n2)
    spec = TheoreticalInstrument.demo_generate_spectrum_line(n2)
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
using Statistics

###############################################################################
# Stacking routines

"""
    combine(cube, [angles]; method=median)

Combine all the frames of a cube using `method`. If `angles` are provided, will use [`derotate!`](@ref) before combining.

Note that with `angles` this is slightly slower than doing `combine(derotate!(cube, angles))`, if you can overwrite the input.

# Examples
```jldoctest
julia> X = ones(2, 3, 3)
2×3×3 Array{Float64,3}:
[:, :, 1] =
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 2] =
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 3] =
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> combine(X)
3×3 Array{Float64,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> combine(X, [0, 90])
1×3×3 Array{Float64,3}:
[:, :, 1] =
 1.0  1.0  1.0

[:, :, 2] =
 NaN  1.0  NaN

[:, :, 3] =
 1.0  1.0  1.0
```
"""
function combine(cube::AbstractArray{T,3}; method = median) where T
    return method(cube, dims = 1)[1, :, :]
end

combine(cube::AbstractArray{T,3}, angles::AbstractVector; method = median) where T = combine(derotate(cube, angles); method = method)

"""
    flatten(cube)

Given a cube of size `(n, x, y)` returns a matrix with size `(n, x * y)` where each row is a flattened image from the cube.

# Examples
```jldoctest
julia> X = ones(3, 2, 2)
3×2×2 Array{Float64,3}:
[:, :, 1] =
 1.0  1.0
 1.0  1.0
 1.0  1.0

[:, :, 2] =
 1.0  1.0
 1.0  1.0
 1.0  1.0

julia> flatten(X)
3×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
```

# See Also
* [`expand`](@ref)
"""
flatten(cube::AbstractArray{T,3}) where T = reshape(cube, size(cube, 1), size(cube, 2) * size(cube, 3))

"""
    expand(matrix)

Given a matrix of size `(n, z)`, returns a cube of size `(n, x, x)` where `x=√z`.

Will throw an error if `z` is not a perfect square.

# Examples
```jldoctest
julia> X = ones(3, 4)
3×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> expand(X)
3×2×2 Array{Float64,3}:
[:, :, 1] =
 1.0  1.0
 1.0  1.0
 1.0  1.0

[:, :, 2] =
 1.0  1.0
 1.0  1.0
 1.0  1.0
```

# See Also
* [`flatten`](@ref)
"""
function expand(mat::AbstractMatrix)
    n, z = size(mat)
    x = sqrt(z)
    isinteger(x) || error("Array of size $((n, x, x)) is not compatible with input matrix of size $(size(mat)).")
    return reshape(mat, n, Int(x), Int(x))
end

###############################################################################
# Rotation routines

using ImageTransformations: imrotate, center

"""
    derotate!(cube, angles)

In-place version of [`derotate`](@ref)
"""
function derotate!(cube::AbstractArray{T,3}, angles::AbstractVector) where {T <: AbstractFloat}
    @inbounds for i in axes(cube, 1)
        frame = @view cube[i, :, :]
        cube[i, :, :] .= imrotate(frame, deg2rad(angles[i]), axes(frame))
    end
    return cube
end

"""
    derotate(cube, angles)

Rotates an array using the given angles in degrees.

This will rotate frame `i` counter-clockwise by the amount `deg2rad(angles[i])`.

# See Also
* [`derotate!`](@ref), [`rotate`](@ref), [`rotate!`](@ref)
"""
function derotate(cube::AbstractArray{T,3}, angles::AbstractVector) where T
    all(angles .≈ 0) && return cube
    return derotate!(deepcopy(cube), angles)
end

"""
    rotate!(cube, angles)

In-place version of [`derotate`](@ref)
"""
function rotate!(cube::AbstractArray{T,3}, angles::AbstractVector) where {T <: AbstractFloat}
    @inbounds for i in axes(cube, 1)
        frame = @view cube[i, :, :]
        cube[i, :, :] .= imrotate(frame, -deg2rad(angles[i]), axes(frame))
    end
    return cube
end

"""
    rotate(cube, angles)

Rotates an array using the given angles in degrees.

This will rotate frame `i` clockwise by the amount `deg2rad(angles[i])`.

# See Also
* [`rotate!`](@ref), [`derotate`](@ref), [`derotate!`](@ref)
"""
function rotate(cube::AbstractArray{T,3}, angles::AbstractVector) where T
    all(angles .≈ 0) && return cube
    return rotate!(deepcopy(cube), angles)
end


text:
let prec = precision(BigFloat), spi = "", digit = 1
    while true
      if digit > lastindex(spi)
        prec *= 2
        setprecision(prec)
        spi = string(big(π))
      end
      print(spi[digit])
      digit += 1
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using LinearAlgebra
using Statistics
using DelimitedFiles;
using jInv.Mesh;
using MAT;
using jInv
using jInv.InverseSolve
using ShapeReconstructionPaLS
using ShapeReconstructionPaLS.PointCloud;
using ParamLevelSet;
using ShapeReconstructionPaLS.ShapeFromSilhouette;
@everywhere using ShapeReconstructionPaLS.Utils;
using Statistics
using Distributed
using SparseArrays
using DelimitedFiles
using Random

ENV["MPLBACKEND"] = "Qt4Agg"

plotting = true;

if plotting
    using PyPlot;
    using jInvVisPyPlot;
    close("all");
end

####
#Parameters to choose if: inversion with dip only/visual hull only/ joint inversion
invertPC = true;
invertVis = false;

#Place new added RBFs according to gradient values or randomly:
locateRBFwithGrads = invertPC;
pad_in_pixels = 25;
n = [100,100,100].+2*pad_in_pixels;
Mesh = getRegularMesh([0.0 1.5 0.0 1.5 0.0 1.5],n);

################################################
### Reading the model and pad ##################
################################################

model = "fandiskpc";
# model = "fandisksmall";
if invertPC
	file_name = string(pwd(),"/../models/",model,".xyz");
	PC_orig = readdlm(file_name);
	Normals = PC_orig[:,4:6];
	PC = PC_orig[:,1:3];
	#PC = PC*Diagonal(0.5./maximum(PC,dims=1)[:]);
	PC = PC .+ 0.5*(Mesh.domain[2] + Mesh.domain[1]);
	##Sort the PC:
	PC  = [ PC Normals];
	PC = PC[sortperm(PC[:, 1]), :];
	writedlm("SortedPC.txt",PC);
	println("PC:",size(PC)," ", PC[end,:]);
	if plotting
		scatter3D(Vector(PC[:,1]),Vector(PC[:,2]),Vector(PC[:,3]),s=0.1,alpha=1.0,color="blue")
	end
end

global m = 0;
if invertVis
	##Read mat file for for visual hull:
	file = matopen(string(pwd(),"/../models/",model,".mat"));
	A = read(file, "B");
	close(file);
	A = convert(Array{Float32,3},A[1:4:end,1:4:end,1:4:end]); # that's 100^3.
	m = zeros(Float64,tuple(n...))
	# m[pad_in_pixels+1:end-pad_in_pixels,pad_in_pixels+1:end-pad_in_pixels,pad_in_pixels+1:end-pad_in_pixels] = A;
	pad = div.(collect(n) .- collect(size(A)),2);
	m[pad[1]+1:end-pad[1],pad[2]+1:end-pad[2]-1,pad[3]+1:end-pad[3]] .= A;
	A=[];
end
# if plotting
    # plotModel(m);
    # savefig("TrueModel.png");
# end

################################################################
### Preparing the synthetic data using mat-free implementation
################################################################
ScreenMesh = getRegularMesh(Mesh.domain[3:6],n[2:3]); ## for vis.
ScreenDomain = ScreenMesh.domain;
n_screen = ScreenMesh.n;

#########################################################################################################
### Point Cloud Data ########################################################################################
#########################################################################################################

npc = 0;
nShots = 0;
noiseTrans = 0.015;
noiseAngles = 5.0;
if invertPC
	println("Data gen resolution:")
	println(n)
	dipDataFilename = string("pointCloudData",model,"_dataRes",n,1);
	trans = [0.0 0 0 ; 0.0 0 0 ]
	theta_phi_dip = [0.0 0.0 ; 0.0 0.0];
	noiseTrans = 0.1;
	noiseAngles = deg2rad(noiseAngles);
	npc = 2;
	prepareSyntheticPointCloudData(PC,Mesh,npc,theta_phi_dip,trans,noiseTrans,noiseAngles,dipDataFilename);
	PC = [];
	#########################################################################################################
	### Read the dipping data ###############################################################################
	#########################################################################################################
	DataPC,P,Normals,trans,theta_phi_dip,npc = readPointCloudDataFile(dipDataFilename);
	### Set up the dip inversion
	P =  Array{Array{Float64,2}}(P);
	Normals = Array{Array{Float64,2}}(Normals);

	if plotting
		close(998);
		figure(998);
		idxs = round(Int64,size(P[1],1)/3)+1:round(Int64,2*size(P[1],1)/3)
		scatter3D(Vector(P[1][idxs,1]),Vector(P[1][idxs,2]),Vector(P[1][idxs,3]),s=0.1,alpha=1.0,color="blue")
		idxs = round(Int64,size(P[2],1)/3)+1:round(Int64,2*size(P[2],1)/3)
		scatter3D(Vector(P[2][idxs,1]),Vector(P[2][idxs,2]),Vector(P[2][idxs,3]),s=0.1,alpha=1.0,color="red")
	end
end

####################################

#################################################################################################################
### Generate and Read the visual data ###############################################################################
#################################################################################################################
noiseAnglesVis = 0.0;
noiseTransVis = 0.0;
samplingBinning = 1;
if invertVis
	nShots = 8;
	noiseSample = 0.0;
	factorShotsCreation = 4;
	theta_phi_vis = deg2rad.(convert(Array{Float64,2}, [rand(0:359,factorShotsCreation*nShots)  rand(0:90,factorShotsCreation*nShots)]));
	VisDataFilename = string("visData_",model,"_dataRes",n,"_noiseSampleAnglesTrans",[noiseSample;noiseAngles;noiseTrans],"_shots",nShots);

	XlocScreen = Mesh.domain[2];
	LocCamera = [-10.0*(Mesh.domain[2] - Mesh.domain[1]) ; 0.0 ;0.0]; ## 0.0 here is the middle of the domain
	noiseCamLoc = 0.0;
	noiseScreenLoc = 0.0;

    noiseAnglesVis = deg2rad(0.0);
	noiseTransVis = 0.0;
	prepareSyntheticVisDataFiles(m,Mesh,ScreenMesh,pad_in_pixels,theta_phi_vis,XlocScreen,LocCamera,
	VisDataFilename,factorShotsCreation*nShots, noiseCamLoc, noiseScreenLoc, noiseAnglesVis,noiseTransVis,MATFree);

	##################################################################################################################
	##################################################################################################################

	n,domain,XScreenLoc,CameraLoc,ScreenDomain,n_screen,DataVis,theta_phi_vis = readVisDataFile(VisDataFilename);
	norms = vec(sum(DataVis,dims = 1));
	idsorted = sortperm(norms);
	theta_phi_vis = theta_phi_vis[idsorted[(end-nShots+1):end],:];
	DataVis = DataVis[:,idsorted[(end-nShots+1):end]];

	ScreenMesh 	= getRegularMesh(ScreenDomain,div.(n_screen,samplingBinning));
	Mesh = getRegularMesh(domain,div.(n,samplingBinning));
	# pad = div(pad,samplingBinning);
	if  invertVis #samplingBinning==2 &&
		#DataVis = coarsenByTwo(DataVis,ScreenMesh.n,nShots);
		if plotting
			figure();
			for kkk = 1:nShots
				dk = DataVis[:,kkk];
				dk = reshape(dk,tuple(ScreenMesh.n...));
					subplot(8,8,kkk);
				plotModel(dk);

			end
			savefig(string("VisTrueModel.png"));
		end
	end
end

function plotVisData(Dc)
if plotting
	close("all");
	close(999);
	figure(999);
	for kkk = 1:nShots
		dk = Dc[:,kkk];
		dk = reshape(dk,tuple(ScreenMesh.n...));
		subplot(2,2,kkk);
		plotModel(dk);
	end
end
end

# function plotPCData(PC)
# if plotting
	# close("all");
	# close(998);
	# figure(998);
	# for kkk = 1:npc
		# dk = Dc[:,kkk];
		# dk = reshape(dk,tuple(ScreenMesh.n...));
		# subplot(2,2,kkk);
		# plotModel(dk);
	# end
# end
# end

samplingBinning = 1;
println("Inv Mesh:")
println(Mesh.n)

method = RBF10Based; methodName="RBF10Based";

n_tup = tuple(Mesh.n...);
nWorkers = nworkers();
misfun = SSDFun; ## least squares

if invertPC
	### Set up the dip inversion
	# 0.0 below is margin : TODO: remove margin.
	pForPC = getPointCloudParam(Mesh,P,Normals,theta_phi_dip,trans,1,method);

	### Create Dip param if we use the GRAD of MATFREE to locate new RBFs
	if(locateRBFwithGrads)
		pForPC_MATFree = getPointCloudParam(Mesh,P,Normals,theta_phi_dip,trans,1,MATFree);
	end

	dobsPC = divideDataToWorkersData(1,DataPC);
	Wd_PC = Array{Array{Float32}}(undef,1);
	for k=1:1
		Wd_PC[k] = ones(Float32,size(dobsPC[k]));
	end
end

if invertVis
	### Set up the vis inversion
	pForVis = getSfSParam(Mesh,theta_phi_vis,zeros(nShots,3),ScreenMesh,XlocScreen,LocCamera,method,nWorkers);
	dobsVis = divideDataToWorkersData(nWorkers,DataVis);
	Wd_Vis = Array{Array{Float32}}(undef,nWorkers);
	for k=1:length(Wd_Vis)
		println("sqrt:",0.1*(Mesh.domain[2]-Mesh.domain[1])*prod(Mesh.h[2:3]));
		Wd_Vis[k] = 1e+1*sqrt((Mesh.domain[2]-Mesh.domain[1])*prod(Mesh.h[2:3]))*ones(Float32,size(dobsVis[k]));
		# Wd_Vis[k] = sqrt(0.3)*ones(Float32,size(dobsVis[k]));

	end
end

if invertVis && invertPC
	n_Moves_all = nShots + npc;
elseif invertPC
	n_Moves_all = npc;
else
	n_Moves_all = nShots;
end

#########################################################################################################
############### Volumetric inversion (pixel-wise) #######################################################
#########################################################################################################
if method == MATBased || method == MATFree
	error("Not relevant for PC.")
else
	nRBF = 2;

	isRBF10 = (method == RBF10BasedSimple1 || method==RBF10BasedSimple2 || method==RBF10Based);
	isSimple = (method==RBFBasedSimple1 || method==RBF10BasedSimple1 || method==RBFBasedSimple2 || method==RBF10BasedSimple2);

	if isRBF10
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
	n_m_simple = numParamOfRBF*nRBF;
	nAll = isSimple ? numParamOfRBF*nRBF : numParamOfRBF*nRBF + n_Moves_all*5;

	modfun = identityMod;

	if isSimple
		modfunForPlotting = (m)->MeshFreeParamLevelSetModelFunc(Mesh,m;computeJacobian = 0,
		sigma = getDefaultHeaviside() ,bf = 1, numParamOfRBF = numParamOfRBF);
	else
		modfunForPlotting = (m)->MeshFreeParamLevelSetModelFunc(Mesh,m[1:(numParamOfRBF*nRBF)];computeJacobian = 0,
		sigma = getDefaultHeaviside() ,bf = 1, numParamOfRBF = numParamOfRBF);
	end

	Iact = 1.0;
	IactPlot = 1.0;
	sback = zero(Float32);
	a = -100000000.0;
	b = 10000000.0;
	boundsHigh = zeros(Float32,nAll) .+ b;
	boundsLow = zeros(Float32,nAll) .+ a;
	# m0 - initial guess. mref - reference for regularization.
	mref = zeros(Float32,nAll);
	m0   = zeros(Float64,n_m_simple);
	mCenters = zeros(Float32,n_tup);
	mref[1:numParamOfRBF:n_m_simple] .= ParamLevelSet.centerHeavySide*2.0;
	if isRBF10
		mref[2:10:n_m_simple] .= 4.0;
		mref[5:10:n_m_simple] .= 4.0;
		mref[7:10:n_m_simple] .= 4.0;
		boundsLow[2:10:n_m_simple] .= 0.05;
		boundsLow[5:10:n_m_simple] .= 0.05;
		boundsLow[7:10:n_m_simple] .= 0.05;
	else
		mref[2:5:n_m_simple] .= 2.0;
		boundsLow[2:5:n_m_simple] .= 0.05;
	end
	mref[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .= (Mesh.domain[1] + Mesh.domain[2])/2;
	mref[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .= (Mesh.domain[3] + Mesh.domain[4])/2;
	mref[numParamOfRBF:numParamOfRBF:n_m_simple] .= (Mesh.domain[5] + Mesh.domain[6])/2;
	m0[:] = mref[1:n_m_simple];
	if !isSimple
		if invertPC
			m0 = wrapRBFparamAndRotationsTranslations(m0,theta_phi_dip,trans);
		end
		if invertVis
			m0 = wrapRBFparamAndRotationsTranslations(m0,theta_phi_vis,zeros(size(theta_phi_vis,1),3));
		end
	end
	mref[:] = m0;
	mref[(n_m_simple+1):(n_m_simple+5*npc)].= 0.0;
	m0[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .+= 0.1*(Mesh.domain[2] - Mesh.domain[1])*randn(nRBF);
	m0[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .+= 0.1*(Mesh.domain[4] - Mesh.domain[3])*randn(nRBF);
	m0[numParamOfRBF:numParamOfRBF:n_m_simple]     .+= 0.1*(Mesh.domain[6] - Mesh.domain[5])*randn(nRBF);

	alpha = 1e+3;

	II = (sparse(1.0I, nAll,nAll));
	II.nzval[(n_m_simple+1):(n_m_simple+5*npc)] .= 1e+4;
	regfun = (m, mref, M)->TikhonovReg(m,mref,M,II);
	# HesPrec = getEmptyRegularizationPreconditioner();
	HesPrec = getSSORCGRegularizationPreconditioner(1.0,1e-2,1);

	if isRBF10
		spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
		regfun 	 		= [regfun;spd_reg];
		alpha 		 	= [alpha;(1e-5)*alpha];
		mref 			= [mref zeros(size(mref))];
	end
end

################################################################################################################
####### Set up inversion #######################################################################################
################################################################################################################

### In case we use MATFREE gradMisfit to locate new RBFs
if(locateRBFwithGrads)
	mask=0;
	Iact_free = 1.0;
	#### USE A BOUND MODEL or bounds
	a_free = 0.0;
	b_free= 1.0;

	pMis_Free = getMisfitParam(pForPC_MATFree, Wd_PC, dobsPC, misfun, Iact_free,sback);
end
if invertVis && invertPC
	if isSimple
		Wd   = [Wd_PC;Wd_Vis];
		dobs = [dobsPC;dobsVis];
		pFor = [pForPC;pForVis];
		pMisRFs = getMisfitParam(pFor, Wd, dobs, misfun, Iact,sback);
	else
		modfun1 = (mAll) -> splitTwoRotationsTranslations(mAll,npc,nShots,1);
		modfun2 = (mAll) -> splitTwoRotationsTranslations(mAll,npc,nShots,2);
		pMisRFs1 = getMisfitParam(pForPC, Wd_PC, dobsPC, misfun, Iact,sback,ones(length(pForPC)),modfun1);
		pMisRFs2 = getMisfitParam(pForVis, Wd_Vis, dobsVis, misfun, Iact,sback,ones(length(pForVis)),modfun2);
		pMisRFs = [pMisRFs1;pMisRFs2];
		pMisRFs1 = 0;
		pMisRFs2 = 0;
	end
	n_Moves_all = nShots + npc;

elseif invertPC
	pMisRFs = getMisfitParam(pForPC, Wd_PC, dobsPC, misfun, Iact,sback);
else
	pMisRFs = getMisfitParam(pForVis, Wd_Vis, dobsVis, misfun, Iact,sback);
end

cgit  = 50; # not used.
maxit = 5;
pcgTol = 1e-2; #not used

@everywhere function updateAngles(pFor::Array{RemoteChannel},theta_phi)
@sync begin
	@async begin
		for k=1:length(pFor)
			pFor[k] = remotecall_fetch(updateAngles,pFor[k].where,pFor[k],theta_phi);
		end
	end
end
return pFor;
end

@everywhere function updateAngles(pFor::RemoteChannel,theta_phi)
	pForMATFREE  = take!(pFor)
	worker_theta_phi = theta_phi[pForMATFREE.workerSubIdxs,:]
	pForMATFREE.theta_phi_rad=copy(worker_theta_phi);
	put!(pFor,pForMATFREE)
	return pFor;
end
@everywhere function updateAngles(pFor::RemoteChannel,theta_phi)
pForMATFREE  = take!(pFor)
worker_theta_phi = theta_phi[pForMATFREE.workerSubIdxs,:]
pForMATFREE.theta_phi_rad=copy(worker_theta_phi);
put!(pFor,pForMATFREE)
return pFor;
end

function myDump(mc::Vector,Dc,iter,pInv,pMis,method="",npc=0,noiseAngle=0,noiseTranslation=0,useVisual=false,nshots=4,plotFun=modfunForPlotting)
	ntup = tuple((pInv.MInv.n)...);
	fullMc = plotFun(mc)[1];
	fullMc = reshape(IactPlot*fullMc[:],ntup);
	# fullMc = reshape(fullMc[:],ntup);
	close(888);
	figure(888);
	#plotModel(fullMc,true,pInv.MInv);
	plotModel(fullMc,includeMeshInfo = true,M_regular = pInv.MInv);
	sleep(1.0);
	savefig(string("It_2nd",iter,"_",ntup,".png"));
	writedlm(string("Params_2nd_iter_",iter,".dat"),convert(Array{Float16},mc));
	writedlm(string("Volume_2nd_iter_",iter,"_",ntup,".dat"),convert(Array{Float16},fullMc));
end

pInv = getInverseParam(Mesh,modfun,regfun,alpha,mref,boundsLow,boundsHigh,
maxStep = 2e-1,pcgMaxIter=cgit,pcgTol=pcgTol,
minUpdate = 0.001, maxIter = maxit,HesPrec=HesPrec);

#### Projected Gauss Newton
mc = m0;
Dc = nothing
pInv.maxIter = 10;
global mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
pInv.maxIter = 20;
myDump(mc,0,1,pInv,0,methodName,npc,noiseAngles,noiseTrans,invertVis,nShots);

if !(method == MATBased || method == MATFree)
	outerIter = 80;
	grad_u = nothing;
	new_RBF_location = [];
	@sync for iterNum=2:outerIter
		new_nRBF = 10;
		println("Inversion using ", nRBF + new_nRBF," basis functions")

		# Dc - the data using our method.
		global Dc = Dc;

		#Lets compute the gradient of the misfit function:
		# locateRBFwithGrads = true;
		print(mc[end-5*n_Moves_all+1:end])
		if(locateRBFwithGrads)
			#(m11,theta_phi1,b1) 	= splitRBFparamAndRotationsTranslations(mc,nRBF,ndips,numParamOfRBF)
			#updateAngles(pForDip_MATFree,theta_phi1);
			pMis_Free = getMisfitParam(pForPC_MATFree, Wd_PC, dobsPC, misfun, Iact,sback);
			sigmaH = getDefaultHeaviside();
			allPoints = [P[1] ; P[2]];
			Dd, = computeMisfit(mc,fetch(pMisRFs[1])); #first element is pmis of PC
			#gradMis = (fetch(Dd[1])-DataPC).^2;
			gradMis = ((Dd)-DataPC).^2;
			sp = sortperm(gradMis,rev = true);
			sp = sp[1:round(Int64,length(sp)/10.0)];
			sp = sp[randperm(length(sp))];
			idxs = sp[1:new_nRBF];
			bwdSubs = [P[1] ; P[2]];
			println("all P size:",size(allPoints));
			Ii = allPoints[idxs,:];
		else
			u = modfunForPlotting(mc)[1];
			#Ii = findall(x -> x >0.2 && x<0.5,u);
			Ii = findall(x -> x >0.2 && x<0.9,u);
			if length(Ii) > new_nRBF
				Ii = Ii[randperm(length(Ii))[1:new_nRBF]];
			end
		end
    #while(Ii in new_RBF_location || Ii.-1 in new_RBF_location || Ii.+1 in new_RBF_location)
    #  Ii = findall(x -> x .>0.2 && x.<0.5,u);
    #  Ii = Ii[randperm(length(Ii))[1:new_nRBF]];
    #end
    	# println("length of I",size(Ii));
		global new_RBF_location = append!(new_RBF_location,Ii);

		if(locateRBFwithGrads)
			# println("Added rbfs at locations:",Ii);
			X1 = Ii[:,1]; X2 = Ii[:,2]; X3 = Ii[:,3];
		else
			println("ii:",Ii)
			println(ind2subv(n_tup,Ii))
			indices = getindex.(ind2subv(n_tup,Ii),[1 2 3 ]);
			println(indices)

			I1 = indices[:,1];
			I2 = indices[:,2];
			I3 = indices[:,3];

			X1 = I1*Mesh.h[1] .- 0.5*Mesh.h[1] .+ Mesh.domain[1];
			X2 = I2*Mesh.h[2] .- 0.5*Mesh.h[2] .+ Mesh.domain[3];
			X3 = I3*Mesh.h[3] .- 0.5*Mesh.h[3] .+ Mesh.domain[5];
		end
		# println("X1:",X1);println("X2:",X2);println("X3:",X3);

		n_mc_old = numParamOfRBF*nRBF;
		n_mc_new = numParamOfRBF*nRBF + size(Ii,1)*numParamOfRBF;
		mc_new = zeros(length(mc) + size(Ii,1)*numParamOfRBF);
		if isRBF10
			mref = pInv.mref[:,1];
		else
			mref = pInv.mref;
		end
		mref_new = zeros(Float32,length(mc) + size(Ii,1)*numParamOfRBF);

		mc_new[1:n_mc_old] = mc[1:n_mc_old];

		mc_new[(n_mc_new+1):end] = mc[(n_mc_old+1):end];
		## Below we set the reference for moves (misalignments) to be 0, assuming they are small
		mref_new[(n_mc_new+1):end] = mref[(n_mc_old+1):end];
		## Below we set the reference for moves (misalignments) to be the previous iteration, assuming they are not so small, and can accumulate
		# mref_new[(n_mc_new+1):end] = mc[(n_mc_old+1):end];

		mc_new[(n_mc_old+1):numParamOfRBF:n_mc_new] .= 0.0; #ParamLevelSet.centerHeavySide - ParamLevelSet.deltaHeavySide;
		mc_new[(n_mc_old+numParamOfRBF-2):numParamOfRBF:n_mc_new] = X1;
		mc_new[(n_mc_old+numParamOfRBF-1):numParamOfRBF:n_mc_new] = X2;
		mc_new[(n_mc_old+numParamOfRBF):numParamOfRBF:n_mc_new] = X3;

		boundsLow = zeros(Float32,length(mc_new)).+a;
		boundsHigh = zeros(Float32,length(mc_new)).+b;
		boundsLow[1:n_mc_old] = pInv.boundsLow[1:n_mc_old];
		boundsHigh[1:n_mc_old] = pInv.boundsHigh[1:n_mc_old];

		eeps =  1.0/(0.01171875*10.0*(0.95^iterNum))
		if isRBF10
			mc_new[(n_mc_old+2):10:n_mc_new] .= eeps #10.0;
			mc_new[(n_mc_old+5):10:n_mc_new] .= eeps#10.0;
			mc_new[(n_mc_old+7):10:n_mc_new] .= eeps#10.0;
			boundsLow[(n_mc_old+2):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+5):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+7):10:n_mc_new] .= 0.05;
		else
			mc_new[(n_mc_old+2):5:n_mc_new] .= 6.0;
			boundsLow[(n_mc_old+2):5:n_mc_new] .= 0.05;
		end
		# boundsLow[end-9:end] .= -0.01;
		# boundsHigh[end-9:end] .= 0.01;
		pInv.boundsLow = boundsLow;
		pInv.boundsHigh = boundsHigh;
		mref_new[1:n_mc_new] = mc_new[1:n_mc_new];
		# mref_new[end-5*n_Moves_all+1:end-5*(n_Moves_all - nShots)].=0.0;
		# the next line tries to make first PC not to move
		mref_new[(n_mc_new+1):(n_mc_new+5)].=0.0;
		global mc = mc_new;
		global nRBF += new_nRBF;

		#II = speye(Float32,length(mc_new));
		len = Int64(length(mc_new));
		IIs = sparse(1.0I,len,len);
		IIs.nzval[(n_mc_new+1):(n_mc_new+5*npc)] .= 1e+4;
		pInv.regularizer = (m, mref, M)->TikhonovReg(m,mref,M,IIs);

		if isRBF10
			spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
			pInv.regularizer  = [pInv.regularizer;spd_reg];
			pInv.mref 			  = [mref_new zeros(length(mref_new))];
			pInv.alpha[1]  = pInv.alpha[1]*0.9;
		else
			pInv.mref = mref_new;
			pInv.alpha[1] *= pInv.alpha[1]*0.9;
		end

		mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
		#mc[end-5*npc+1:end] .= 0.0;
		myDump(mc,0,iterNum,pInv,0,methodName,0,noiseAngles,noiseTrans,invertVis);

		if iterNum==outerIter
			pInv.alpha *= 0.1;
			pInv.maxIter = 200;
			mc, = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
			myDump(mc,0,iterNum+1,pInv,0,methodName,0,noiseAngles,noiseTrans,invertVis);
		end

	end
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
#####################################################################
# BIM
#=
Building Information Modeling is more than just shapes.
Each BIM element is connected to other BIM elements. A wall is connected to its
windows and doors. A floor is connected to its walls, and stairs. Stairs connect
different floors, etc.

This graph of objects needs to be build programmatically and that is one problem.
The other problem is portability, as we want this to work in backends that are
not BIM tools.

A third problem is related to level of detail. Even in the same backend, we might
need to represent BIM elements at different levels of detail.

We will start with the simplest of BIM elements, namely, levels, slabs, walls,
windows and doors.
=#

abstract type BIMElement <: Proxy end
abstract type Measure <: BIMElement end
BIMElements = Vector{<:BIMElement}

# Level
@defproxy(level, Measure, height::Real=0, elements::BIMElements=BIMElement[])
levels_cache = Dict{Real,Level}()
maybe_replace(level::Level) = get!(levels_cache, level.height, level)

current_levels() = values(level_cache)
default_level = Parameter{Level}(level())
default_level_to_level_height = Parameter{Real}(3)
upper_level(lvl=default_level(), height=default_level_to_level_height()) = level(lvl.height + height, backend=backend(lvl))
Base.:(==)(l1::Level, l2::Level) = l1.height == l2.height

#default implementation
realize(b::Backend, s::Level) = s.height

export all_levels, default_level, default_level_to_level_height, upper_level

#=
@defproxy(polygonal_mass, Shape3D, points::Locs, height::Real)
@defproxy(rectangular_mass, Shape3D, center::Loc, width::Real, len::Real, height::Real)

@defproxy(column, Shape3D, center::Loc, bottom_level::Any, top_level::Any, family::Any)
=#

#=

We need to provide defaults for a lot of things. For example, we want to specify
a wall that goes through a path without having to specify the kind of wall or its
thickness and height.

This means that, apart from the wall's path, all other wall features will come
from defaults. The base height will be determined by the current level and the
wall height by the current level-to-level height. Finally, the wall thickness,
constituent parts, thermal characteristics, and so on will come from the wall
defaults.  In the case of wall, we will assume that current_wall_defaults() is
a parameter that contains a set of wall parameters.  As an example of use, we
might have:

current_wall_defaults(wall_defaults(thickness=10))

Another option is the definition of different defaults:

thick_wall_defaults = wall_defaults(thickness=10)
thin_wall_defaults = wall_defaults(thickness=5)

which then can be make current:

current_wall_defaults(thin_wall_defaults)

In most cases, the defaults are not just one value, but a bunch of them. For a
beam, we might have:

standard_beam = beam_defaults(width=10, height=20)
current_beam_defaults(standard_beam)

Another useful feature is the ability to adapt defaults. For example:

current_beam_defaults(beam_with(standard_beam, height=20))

Finally, defaults can be created for anything. For example, in a building, we
might want to define a bunch of parameters that are relevant. The syntax is as
follows:

@defaults(building,
    width::Real=20,
    length::Real=30,
    height::Real=50)

In order to access these defaults, we can use the following:

current_building_defaults().width

In some cases, defaults are supported by the backend itself. For example, in
Revit, a wall can be specified using a family. In order to realize the wall
defaults in the current backend, we need to map from the wall parameters to the
corresponding BIM family parameters. This mapping must be described in a
different structure.

For example, a beam element might have a section with a given width and height
but, in Revit, a beam element such as "...\\Structural Framing\\Wood\\M_Timber.rfa"
has, as parameters, the dimensions b and d.  This means that we need a map, such
as Dict(:width => "b", :height => "d")))). So, for a Revit family, we might use:

RevitFamily(
    "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
    Dict(:width => "b", :height => "d"))

To make things more interesting, some families might require instantiation on
different levels. For example, a circular window family in Revit needs to be
loaded using RVTLoadFamily, which requires the name of the family, then it needs
to be instantiated with RVTFamilyElement, which requires the radius of the window,
and finally needs to be inserted on the wall, using RVTInsertWindow, which requires
the opening angle. This might be different for different families, so we need a
flexible way of using the parameters. One hipotesis is to specify those different
moments as different dictionaries.

RevitFamily(
  "C:\\ProgramData\\Autodesk\\RVT 2019\\Libraries\\US Metric\\Windows\\CIRCULAR WINDOW.rfa",
  ("radius_window" => :radius),
  ("angle_window" => :opening_angle))

If no parameters are needed on a particular phases, we might use an empty dictionary
to describe that phase. Given the typical Revit families, defaults parameters seem
to be the best approach here. This means that the previous beam family might be
equivalent to

RevitFamily(
  "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
  ("b" => :width, "d" => :height),
  ())

However, the same beam might have a different mapping in a different backend.
This means that we need another mapping to support different backends. One
possibility is to use something similar to:

backend_family(
    revit => RevitFamily(
        "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
        ("b" => :width, "d" => :height)),
    archicad => ArchiCADFamily(
        "BeamElement",
        ("size_x" => :width, "size_y" => :height)),
    autocad => AutoCADFamily())

Then, we need an operation that instantiates a family. This can be done on two different
levels: (1) from a backend-specific family (e.g., RevitFamily), for example:

beam_family = RevitFamily(
    "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
    ("b" => :width, "d" => :height))

current_beam_defaults(beam_family_instance(beam_family, width=10, height=20)

or from a generic backend family, for example:

beam_family = backend_family(
    revit => RevitFamily(
        "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
        ("b" => :width, "d" => :height)),
    archicad => ArchiCADFamily(
        "BeamElement"
        ("size_x" => :width, "size_y" => :height)),
    autocad => AutoCADFamily())

current_beam_defaults(beam_family_instance(beam_family, width=10, height=20)

In this last case, the generic family will use the current_backend value to identify
which family to use.

Another important feature is the use of a delegation-based implementation for
family instances. This means that we might do

current_beam_defaults(beam_family_instance(current_beam_defaults(), width=20)

to instantiate a family that uses, by default, the same parameter values used by
another family instance.

=#

abstract type Family <: Proxy end
abstract type FamilyInstance <: Family end

family(f::Family) = f
family(f::FamilyInstance) = f.family

#HACK Using Dict instead of IdDict just because get is not fully implemented for IdDict
#FIXME after updating to Julia 1.2
macro deffamily(name, parent, fields...)
  name_str = string(name)
  abstract_name = esc(Symbol(string))
  struct_name = esc(Symbol(string(map(uppercasefirst,split(name_str,'_'))...)))
  field_names = map(field -> field.args[1].args[1], fields)
  field_types = map(field -> field.args[1].args[2], fields)
  field_inits = map(field -> field.args[2], fields)
  field_renames = map(esc ∘ Symbol ∘ uppercasefirst ∘ string, field_names)
  field_replacements = Dict(zip(field_names, field_renames))
  struct_fields = map((name,typ) -> :($(name) :: $(typ)), field_names, field_types)
#  opt_params = map((name,typ,init) -> :($(name) :: $(typ) = $(init)), field_renames, field_types, field_inits)
#  key_params = map((name,typ,rename) -> :($(name) :: $(typ) = $(rename)), field_names, field_types, field_renames)
#  mk_param(name,typ) = Expr(:kw, Expr(:(::), name, typ))
  mk_param(name,typ,init) = Expr(:kw, Expr(:(::), name, typ), init)
  opt_params = map(mk_param, field_renames, field_types, map(init -> replace_in(init, field_replacements), field_inits))
  key_params = map(mk_param, field_names, field_types, field_renames)
  instance_params = map(mk_param, field_names, field_types, map(name -> :(family.$(name)), field_names))
  constructor_name = esc(name)
  instance_name = esc(Symbol(name_str, "_element")) #"_instance")) beam_family_element or beam_family_instance?
  default_name = esc(Symbol("default_", name_str))
  predicate_name = esc(Symbol("is_", name_str))
  selector_names = map(field_name -> esc(Symbol(name_str, "_", string(field_name))), field_names)
  quote
    export $(constructor_name), $(instance_name), $(default_name), $(predicate_name), $(struct_name)
    struct $struct_name <: $parent
      $(struct_fields...)
      based_on::Union{Family, Nothing}
      implemented_as::IdDict{<:Backend, <:Family}
      ref::Parameter{Any}
    end
    $(constructor_name)($(opt_params...);
                        $(key_params...),
                        based_on=nothing,
                        implemented_as=IdDict{Backend, Family}()) =
      $(struct_name)($(field_names...), based_on, implemented_as, Parameter{Any}(nothing))
    $(instance_name)(family:: Family, implemented_as=copy(family.implemented_as); $(instance_params...)) =
      $(struct_name)($(field_names...), family, implemented_as, Parameter{Any}(nothing))
    $(default_name) = Parameter($(constructor_name)())
    $(predicate_name)(v::$(struct_name)) = true
    $(predicate_name)(v::Any) = false
#    $(map((selector_name, field_name) -> :($(selector_name)(v::$(struct_name)) = v.$(field_name)),
#          selector_names, field_names)...)
    Khepri.meta_program(v::$(struct_name)) =
        Expr(:call, $(Expr(:quote, name)), $(map(field_name -> :(meta_program(v.$(field_name))), field_names)...))
  end
end

# When dispatching a BIM operation to a backend, we also need to dispatch the family

backend_family(b::Backend, family::Family) =
# replace this with next fragment after updating Julia
  if haskey(family.implemented_as, b)
    family.implemented_as[b]
  else
    family.based_on == nothing ? # this is not a family_element (nor a derivation of a family_element)
      error("Family $(family) is missing the implementation for backend $(b)") :
      backend_family(b, family.based_on)
  end
#=  get(family.implemented_as, b) do
    family.based_on == nothing ? # this is not a family_element (nor a derivation of a family_element)
      error("Family $(family) is missing the implementation for backend $(b)") :
      backend_family(b, family.based_on)
  end
=#

copy_struct(s::T) where T = T([getfield(s, k) for k ∈ fieldnames(T)]...)

# Backends will install their own families on top of the default families, e.g.,
# set_backend_family(default_beam_family(), revit, revit_beam_family)
set_backend_family(family::Family, backend::Backend, backend_family::Family) =
  begin
    family.implemented_as[backend]=backend_family
    family.ref(nothing) #force recreation
  end

# Finally, we can implement a generic backend caching mechanism for families

realize(b::Backend, f::Family) =
  if f.ref() == nothing
    f.ref(backend_get_family_ref(b, f, backend_family(b, f)))
  else
    f.ref()
  end

export backend_family, set_backend_family

@deffamily(slab_family, Family,
    thickness::Real=0.2,
    coating_thickness::Real=0.0)

slab_family_elevation(b::Backend, family::SlabFamily) =
  family.coating_thickness - family.thickness
slab_family_thickness(b::Backend, family::SlabFamily) =
  family.coating_thickness + family.thickness

@defproxy(slab, Shape3D, contour::ClosedPath=rectangular_path(),
          level::Level=default_level(), family::SlabFamily=default_slab_family(),
          openings::Vector{<:ClosedPath}=ClosedPath[])

# Default implementation: dispatch on the slab elements
realize(b::Backend, s::Slab) =
    realize_slab(b, s.contour, s.openings, s.level, s.family)

realize_slab(b::Backend, contour::ClosedPath, holes::Vector{<:ClosedPath}, level::Level, family::Family) =
    let base = vz(level.height + slab_family_elevation(b, family)),
        thickness = slab_family_thickness(b, family)
        # Change this to a better named protocol?
        backend_slab(b, translate(contour, base), map(c -> translate(c, base), holes), thickness, family)
    end
#
export add_slab_opening
add_slab_opening(s::Slab=required(), contour::ClosedPath=circular_path()) =
    let b = backend(s)
        push!(s.openings, contour)
        if realized(s)
            set_ref!(s, realize_slab_openings(b, s, ref(s), [contour]))
        end
        s
    end

realize_slab_openings(b::Backend, s::Slab, s_ref, openings) =
    let s_base_height = s.level.height,
        s_thickness = slab_family_thickness(b, s.family)
        for opening in openings
            op_path = translate(opening, vz(s_base_height-1.1*s_thickness))
            op_ref = ensure_ref(b, backend_slab(b, op_path, s_thickness*1.2))
            s_ref = ensure_ref(b, subtract_ref(b, s_ref, op_ref))
        end
        s_ref
    end

# Roof

@deffamily(roof_family, Family,
    thickness::Real=0.2,
    coating_thickness::Real=0.0)

slab_family_elevation(b::Backend, family::RoofFamily) = 0
slab_family_thickness(b::Backend, family::RoofFamily) =
  family.coating_thickness + family.thickness

@defproxy(roof, Shape3D, contour::ClosedPath=rectangular_path(),
          level::Level=default_level(), family::RoofFamily=default_roof_family(),
          openings::Vector{<:ClosedPath}=ClosedPath[])

realize(b::Backend, s::Roof) =
    realize_slab(b, s.contour, s.openings, s.level, s.family)

# Panel

@deffamily(panel_family, Family,
    thickness::Real=0.02)

@defproxy(panel, Shape3D, vertices::Locs=Loc[], level::Any=default_level(), family::Any=default_panel_family())

#TODO Pass the provided backend
realize(b::Backend, s::Panel) =
  let #p1 = s.vertices[1],
      #p2 = s.vertices[2],
      #p3 = s.vertices[3],
      #n = vz(s.family.thickness, cs_from_o_vx_vy(p1, p2-p1, p3-p1))
      verts = in_world.(s.vertices)
      n = vertices_normal(verts)*(s.family.thickness/2)
    ref(irregular_prism(map(p -> p - n, verts), n*2))
  end

#=

A wall contains doors and windows

=#

# Wall

@deffamily(wall_family, Family,
    thickness::Real=0.2)

@defproxy(wall, Shape3D, path::Path=rectangular_path(),
          bottom_level::Level=default_level(),
          top_level::Level=upper_level(convert(Level, bottom_level)),
          family::WallFamily=default_wall_family(),
          offset::Real=is_closed_path(path) ? 1/2 : 0, # offset is relative to the thickness
          doors::Shapes=Shape[], windows::Shapes=Shape[])
#=
Walls can be joined. That is very important because the wall needs to have
uniform thickness along the entire path.
=#
export join_walls
join_walls(wall1, wall2) =
  if wall1.bottom_level != wall2.bottom_level
    error("Walls with different bottom levels")
  elseif wall1.top_level != wall2.top_level
    error("Walls with different top levels")
  elseif wall1.family != wall2.family
    error("Walls with different families")
  elseif wall1.offset != wall2.offset
    error("Walls with different offsets")
  else
    let w = wall(join_paths(wall1.path, wall2.path),
                 wall1.bottom_level, wall1.top_level,
                 wall1.family, wall1.offset),
        len = path_length(wall1.path)
      for (es,l) in ((wall1.doors, 0), (wall2.doors, len))
        for e in es
          add_door(w, e.loc+vx(l), e.family)
        end
      end
      for (es,l) in ((wall1.windows, 0), (wall2.windows, len))
        for e in es
          add_window(w, e.loc+vx(l), e.family)
        end
      end
      for w in (wall1, wall2)
        delete_shapes(w.doors)
        delete_shapes(w.windows)
        delete_shape(w)
      end
      w
    end
  end

join_walls(walls...) =
  reduce(join_walls, walls)

# Right and Left considering observer looking along with curve direction
# a non-closed wall should have a wall offset of zero and a r_thickness and a l_thickness of 1/2*thickness
# a closed wall should have a wall offset of 1/2 and a r_thickness of zero and a l_thickness of 1*thickness
# Thus, we have:
r_thickness(offset, thickness) = (1/2 - offset)*thickness
l_thickness(offset, thickness) = (1/2 + offset)*thickness

r_thickness(w::Wall) = r_thickness(w.offset, w.family.thickness)
l_thickness(w::Wall) = l_thickness(w.offset, w.family.thickness)

# Door

@deffamily(door_family, Family,
  width::Real=1.0,
  height::Real=2.0,
  thickness::Real=0.05)

@defproxy(door, Shape3D, wall::Wall=required(), loc::Loc=u0(), flip_x::Bool=false, flip_y::Bool=false, family::DoorFamily=default_door_family())

# Window

@deffamily(window_family, Family,
  width::Real=1.0,
  height::Real=2.0,
  thickness::Real=0.05)

@defproxy(window, Shape3D, wall::Wall=required(), loc::Loc=u0(), flip_x::Bool=false, flip_y::Bool=false, family::WindowFamily=default_window_family())

# Default implementation
realize(b::Backend, w::Wall) =
  realize_wall_openings(b, w, realize_wall_no_openings(b, w), [w.doors..., w.windows...])

realize_wall_no_openings(b::Backend, w::Wall) =
  let w_base_height = w.bottom_level.height,
      w_height = w.top_level.height - w_base_height,
      w_path = translate(w.path, vz(w_base_height)),
      r_thickness = r_thickness(w),
      l_thickness = l_thickness(w)
    ensure_ref(b, backend_wall(b, w_path, w_height, l_thickness, r_thickness, w.family))
  end

realize_wall_openings(b::Backend, w::Wall, w_ref, openings) =
  let w_base_height = w.bottom_level.height,
      w_height = w.top_level.height - w_base_height,
      w_path = translate(w.path, vz(w_base_height)),
      r_thickness = r_thickness(w),
      l_thickness = l_thickness(w)
    for opening in openings
      w_ref = realize_wall_opening(b, w_ref, w_path, l_thickness, r_thickness, opening, w.family)
      ref(opening)
    end
    w_ref
  end

realize_wall_opening(b::Backend, w_ref, w_path, l_thickness, r_thickness, op, family) =
  let op_base_height = op.loc.y,
      op_height = op.family.height,
      op_path = translate(subpath(w_path, op.loc.x, op.loc.x + op.family.width), vz(op_base_height)),
      op_ref = ensure_ref(b, backend_wall(b, op_path, op_height, l_thickness, r_thickness, family))
    ensure_ref(b, subtract_ref(b, w_ref, op_ref))
  end

realize(b::Backend, s::Union{Door, Window}) =
  let base_height = s.wall.bottom_level.height + s.loc.y,
      height = s.family.height,
      subpath = translate(subpath(s.wall.path, s.loc.x, s.loc.x + s.family.width), vz(base_height)),
      r_thickness = r_thickness(s.wall),
      l_thickness = l_thickness(s.wall),
      thickness = s.family.thickness
    backend_wall(b, subpath, height, (l_thickness - r_thickness + thickness)/2, (r_thickness - l_thickness + thickness)/2, s.family)
  end

##

export add_door
add_door(w::Wall=required(), loc::Loc=u0(), family::DoorFamily=default_door_family()) =
  backend_add_door(backend(w), w, loc, family)

backend_add_door(b::Backend, w::Wall, loc::Loc, family::DoorFamily) =
    let d = door(w, loc, family=family)
        push!(w.doors, d)
        if realized(w)
            set_ref!(w, realize_wall_openings(b, w, ref(w), [d]))
        end
        w
    end

#
export add_window
add_window(w::Wall=required(), loc::Loc=u0(), family::WindowFamily=default_window_family()) =
  backend_add_window(backend(w), w, loc, family)

backend_add_window(b::Backend, w::Wall, loc::Loc, family::WindowFamily) =
  let d = window(w, loc, family=family)
    push!(w.windows, d)
    if realized(w)
        set_ref!(w, realize_wall_openings(b, w, ref(w), [d]))
    end
    w
  end

#=
A curtain wall is a special kind of wall that is made of a frame with windows.
=#

@deffamily(curtain_wall_frame_family, Family,
  width::Real=0.1,
  depth::Real=0.1,
  depth_offset::Real=0.25)

@deffamily(curtain_wall_family, Family,
  max_panel_dx::Real=1,
  max_panel_dy::Real=2,
  panel::PanelFamily=panel_family(thickness=0.05),
  boundary_frame::CurtainWallFrameFamily=
    curtain_wall_frame_family(width=0.1,depth=0.1,depth_offset=0.25),
  mullion_frame::CurtainWallFrameFamily=
    curtain_wall_frame_family(width=0.08,depth=0.09,depth_offset=0.2),
  transom_frame::CurtainWallFrameFamily=
    curtain_wall_frame_family(width=0.06,depth=0.1,depth_offset=0.11))

@defproxy(curtain_wall, Shape3D,
          path::Path=rectangular_path(),
          bottom_level::Level=default_level(),
          top_level::Level=upper_level(bottom_level),
          family::CurtainWallFamily=default_curtain_wall_family(),
          offset::Real=0.0)
curtain_wall(p0::Loc, p1::Loc;
     bottom_level::Level=default_level(),
     top_level::Level=upper_level(bottom_level),
     family::CurtainWallFamily=default_curtain_wall_family(),
     offset::Real=0.0) =
  curtain_wall([p0, p1], bottom_level=bottom_level, top_level=top_level,
         family=family, offset=offset)

realize(b::Backend, s::CurtainWall) =
  let th = s.family.panel.thickness,
      bfw = s.family.boundary_frame.width,
      bfd = s.family.boundary_frame.depth,
      bfdo = s.family.boundary_frame.depth_offset,
      mfw = s.family.mullion_frame.width,
      mfd = s.family.mullion_frame.depth,
      mdfo = s.family.mullion_frame.depth_offset,
      tfw = s.family.transom_frame.width,
      tfd = s.family.transom_frame.depth,
      tfdo = s.family.transom_frame.depth_offset,
      path = curtain_wall_path(b, s, s.family.panel),
      path_length = path_length(path),
      bottom = level_height(s.bottom_level),
      top = level_height(s.top_level),
      height = top - bottom,
      x_panels = ceil(Int, path_length/s.family.max_panel_dx),
      y_panels = ceil(Int, height/s.family.max_panel_dy),
      refs = []
    push!(refs, backend_curtain_wall(b, s, subpath(path, bfw, path_length-bfw), bottom+bfw, height-2*bfw, th, :panel))
    push!(refs, backend_curtain_wall(b, s, path, bottom, bfw, bfd, :boundary_frame))
    push!(refs, backend_curtain_wall(b, s, path, top-bfw, bfw, bfd, :boundary_frame))
    push!(refs, backend_curtain_wall(b, s, subpath(path, 0, bfw), bottom+bfw, height-2*bfw, bfd, :boundary_frame))
    push!(refs, backend_curtain_wall(b, s, subpath(path, path_length-bfw, path_length), bottom+bfw, height-2*bfw, bfd, :boundary_frame))
    for i in 1:y_panels-1
      l = height/y_panels*i
      sub = subpath(path, bfw, path_length-bfw)
      push!(refs, backend_curtain_wall(b, s, sub, bottom+l-tfw/2, tfw, tfd, :transom_frame))
    end
    for i in 1:x_panels-1
      l = path_length/x_panels*i
      push!(refs, backend_curtain_wall(b, s, subpath(path, l-mfw/2, l+mfw/2), bottom+bfw, height-2*bfw, mfd, :mullion_frame))
    end
    [ensure_ref(b,r) for r in refs]
  end

# By default, curtain wall panels are planar
curtain_wall_path(b::Backend, s::CurtainWall, panel_family::Family) =
  s.path
curtain_wall_path(b::Backend, s::CurtainWall, panel_family::PanelFamily) =
  let path_length = path_length(s.path),
      x_panels = ceil(Int, path_length/s.family.max_panel_dx),
      pts = map(t->in_world(location_at_length(s.path, t)),
                division(0, path_length, x_panels))
    polygonal_path(pts)
  end

backend_curtain_wall(b::Backend, s, path::Path, bottom::Real, height::Real, thickness::Real, kind::Symbol) =
  backend_wall(b, translate(path, vz(bottom)), height, thickness/2, thickness/2, getproperty(s.family, kind))

#
# We need to redefine the default method (maybe add an option to the macro to avoid defining the meta_program)
# This needs to be fixed for windows
#=
meta_program(w::Wall) =
    if isempty(w.doors)
        Expr(:call, :wall,
             meta_program(w.path),
             meta_program(w.bottom_level),
             meta_program(w.top_level),
             meta_program(w.family))
    else
        let door = w.doors[1]
            Expr(:call, :add_door,
                 meta_program(wall(w.path, w.bottom_level, w.top_level, w.family, w.doors[2:end], w.windows)),
                 meta_program(door.loc),
                 meta_program(door.family))
        end
    end
=#

# Beam
# Beams are mainly horizontal elements. By default, a beam is aligned along its top axis
@deffamily(beam_family, Family,
#  width::Real=1.0,
#  height::Real=2.0,
  profile::ClosedPath=top_aligned_rectangular_profile(1, 2))
#beam_family(Width::Real=1.0, Height::Real=2.0; width=Width, height=Height) =
#  beam_family(rectangular_path(xy(-width/2,-height), width, height))

@defproxy(beam, Shape3D, cb::Loc=u0(), h::Real=1, angle::Real=0, family::BeamFamily=default_beam_family())
beam(cb::Loc, ct::Loc, Angle::Real=0, Family::BeamFamily=default_beam_family(); angle::Real=Angle, family::BeamFamily=Family) =
    let (c, h) = position_and_height(cb, ct)
      beam(c, h, angle, family)
    end

# Column
# Columns are mainly vertical elements. A column has its center axis aligned with a line defined by two points

@deffamily(column_family, Family,
    #width::Real=1.0,
    #height::Real=2.0,
  profile::ClosedPath=rectangular_profile(0.2, 0.2))

@defproxy(free_column, Shape3D, cb::Loc=u0(), h::Real=1, angle::Real=0, family::ColumnFamily=default_column_family())
free_column(cb::Loc, ct::Loc, Angle::Real=0, Family::ColumnFamily=default_column_family(); angle::Real=Angle, family::ColumnFamily=Family) =
    let (c, h) = position_and_height(cb, ct)
      free_column(c, h, angle, family)
    end

@defproxy(column, Shape3D, cb::Loc=u0(), angle::Real=0,
  bottom_level::Level=default_level(), top_level::Level=upper_level(bottom_level),
  family::ColumnFamily=default_column_family())

# Default implementation dispatches to the section profile
realize(b::Backend, s::Beam) =
  realize_beam_profile(b, s, s.family.profile, s.cb, s.h)

realize(b::Backend, s::FreeColumn) =
  realize_beam_profile(b, s, s.family.profile, s.cb, s.h)

realize(b::Backend, s::Column) =
  let base_height = s.bottom_level.height,
      height = s.top_level.height - base_height
    realize_beam_profile(b, s, s.family.profile, add_z(s.cb, base_height), height)
  end

# Tables and chairs

@deffamily(table_family, Family,
    length::Real=1.6,
    width::Real=0.9,
    height::Real=0.75,
    top_thickness::Real=0.05,
    leg_thickness::Real=0.05)

@deffamily(chair_family, Family,
    length::Real=0.4,
    width::Real=0.4,
    height::Real=1.0,
    seat_height::Real=0.5,
    thickness::Real=0.05)

@deffamily(table_chair_family, Family,
    table_family::TableFamily=default_table_family(),
    chair_family::ChairFamily=default_chair_family(),
    chairs_top::Int=1,
    chairs_bottom::Int=1,
    chairs_right::Int=2,
    chairs_left::Int=2,
    spacing::Real=0.7)

@defproxy(table, Shape3D, loc::Loc=u0(), angle::Real=0, level::Level=default_level(), family::TableFamily=default_table_family())

realize(b::Backend, s::Table) =
    backend_rectangular_table(b, add_z(s.loc, s.level.height), s.angle, s.family)

@defproxy(chair, Shape3D, loc::Loc=u0(), angle::Real=0, level::Level=default_level(), family::ChairFamily=default_chair_family())

realize(b::Backend, s::Chair) =
    backend_chair(b, add_z(s.loc, s.level.height), s.angle, s.family)

@defproxy(table_and_chairs, Shape3D, loc::Loc=u0(), angle::Real=0, level::Level=default_level(), family::TableChairFamily=default_table_chair_family())

realize(b::Backend, s::TableAndChairs) =
    backend_rectangular_table_and_chairs(b, add_z(s.loc, s.level.height), s.angle, s.family)

# Lights

@defproxy(pointlight, Shape3D, loc::Loc=z(3), color::RGB=rgb(255,255,255), range::Real=10, intensity::Real=4, level::Level=default_level())

realize(b::Backend, s::Pointlight) =
    backend_pointlight(b, add_z(s.loc, s.level.height), s.color, s.range, s.intensity)

@defproxy(spotlight, Shape3D, loc::Loc=z(3), dir::Vec=vz(-1), hotspot::Real=pi/4, falloff::Real=pi/3)

realize(b::Backend, s::Spotlight) =
    backend_spotlight(b, s.loc, s.dir, s.hotspot, s.falloff)

@defproxy(ieslight, Shape3D, file::String=required(), loc::Loc=z(3), dir::Vec=vz(-1), alpha::Real=0, beta::Real=0, gamma::Real=0)

realize(b::Backend, s::Ieslight) =
    backend_ieslight(b, s.file, s.loc, s.dir, s.alpha, s.beta, s.gamma)

#################################

@deffamily(truss_node_family, Family,
    radius::Real=0.2,
    support::Any=false) #(Option node_support)

@deffamily(truss_bar_family, Family,
    radius::Real=0.03)

@defproxy(truss_node, Shape3D, p::Loc=u0(), family::TrussNodeFamily=default_truss_node_family())
@defproxy(truss_bar, Shape3D, p0::Loc=u0(), p1::Loc=u0(), angle::Real=0, family::TrussBarFamily=default_truss_bar_family())

realize(b::Backend, s::TrussNode) = sphere(s.p, s.family.radius)
realize(b::Backend, s::TrussBar) = cylinder(s.p0, s.family.radius, s.p1)

###################################
# BIM
@defop all_levels()
@defop all_walls()
@defop all_walls_at_level(level)