text:
const MOI  = Gurobi.MOI
const MOIT = MOI.Test

const GUROBI_ENV = Gurobi.Env()
const OPTIMIZER = MOI.Bridges.full_bridge_optimizer(
    # Note: we set `DualReductions = 0` so that we never return
    # `INFEASIBLE_OR_UNBOUNDED`.
    Gurobi.Optimizer(
        GUROBI_ENV,
        OutputFlag = 0,
        DualReductions = 0,
        QCPDual = 1,
        InfUnbdInfo = 1,
    ),
    Float64
)

const CONFIG = MOIT.TestConfig()

@testset "Unit Tests" begin
    MOIT.basic_constraint_tests(OPTIMIZER, CONFIG; exclude = [
        (MOI.VectorOfVariables, MOI.SecondOrderCone),
        (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone),
        (MOI.VectorOfVariables, MOI.GeometricMeanCone),
        (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone),
        (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone),
        (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone),
        (MOI.VectorQuadraticFunction{Float64}, MOI.SecondOrderCone),
        (MOI.VectorQuadraticFunction{Float64}, MOI.RotatedSecondOrderCone),
        (MOI.VectorQuadraticFunction{Float64}, MOI.GeometricMeanCone),
    ])
    # TODO(odow): bugs deleting SOC variables. See also the
    # `delete_soc_variables` test.
    MOIT.basic_constraint_tests(
        OPTIMIZER,
        CONFIG;
        include = [
            (MOI.VectorOfVariables, MOI.SecondOrderCone),
            (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone),
            (MOI.VectorOfVariables, MOI.GeometricMeanCone),
            (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone),
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone),
            (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone),
            (MOI.VectorQuadraticFunction{Float64}, MOI.SecondOrderCone),
            (MOI.VectorQuadraticFunction{Float64}, MOI.RotatedSecondOrderCone),
            (MOI.VectorQuadraticFunction{Float64}, MOI.GeometricMeanCone),
        ],
        delete = false
    )
    MOIT.unittest(OPTIMIZER, MOIT.TestConfig(atol=1e-6), [
        # TODO(odow): bug! We can't delete a vector of variables if one is in
        # a second order cone.
        "delete_soc_variables",
        # TODO(odow): implement number of threads.
        "number_threads",
    ])
    MOIT.modificationtest(OPTIMIZER, CONFIG)
end

@testset "Linear tests" begin
    MOIT.contlineartest(OPTIMIZER, MOIT.TestConfig(basis = true), [
        # This requires an infeasiblity certificate for a variable bound.
        "linear12"
    ])
    MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false))
end

@testset "Quadratic tests" begin
    MOIT.contquadratictest(OPTIMIZER, MOIT.TestConfig(atol=1e-3, rtol=1e-3), [
        "ncqcp"  # Gurobi doesn't support non-convex problems.
    ])
end

@testset "Conic tests" begin
    MOIT.lintest(OPTIMIZER, CONFIG)
    MOIT.soctest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3), ["soc3"])
    MOIT.soc3test(
        OPTIMIZER,
        MOIT.TestConfig(duals = false, infeas_certificates = false, atol = 1e-3)
    )
    MOIT.rsoctest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=5e-3))
    MOIT.geomeantest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3))
end

@testset "Integer Linear tests" begin
    MOIT.intlineartest(OPTIMIZER, CONFIG, [
        # Indicator sets not supported.
        "indicator1", "indicator2", "indicator3", "indicator4"
    ])
end

@testset "ModelLike tests" begin

    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "Gurobi"

    @testset "default_objective_test" begin
        MOIT.default_objective_test(OPTIMIZER)
    end

    @testset "default_status_test" begin
        MOIT.default_status_test(OPTIMIZER)
    end

    @testset "nametest" begin
        MOIT.nametest(OPTIMIZER)
    end

    @testset "validtest" begin
        MOIT.validtest(OPTIMIZER)
    end

    @testset "emptytest" begin
        MOIT.emptytest(OPTIMIZER)
    end

    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(OPTIMIZER)
    end

    @testset "copytest" begin
        MOIT.copytest(
            OPTIMIZER,
            MOI.Bridges.full_bridge_optimizer(Gurobi.Optimizer(GUROBI_ENV), Float64)
        )
    end

    @testset "scalar_function_constant_not_zero" begin
        MOIT.scalar_function_constant_not_zero(OPTIMIZER)
    end

    @testset "start_values_test" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag = 0)
        x = MOI.add_variables(model, 2)
        @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimalStart(), x[1]) === nothing
        @test MOI.get(model, MOI.VariablePrimalStart(), x[2]) === nothing
        Gurobi._update_if_necessary(model)
        @test Gurobi.get_dblattrelement(
            model.inner, "Start", Gurobi._info(model, x[1]).column
        ) == Gurobi.GRB_UNDEFINED
        MOI.set(model, MOI.VariablePrimalStart(), x[1], 1.0)
        MOI.set(model, MOI.VariablePrimalStart(), x[2], nothing)
        @test MOI.get(model, MOI.VariablePrimalStart(), x[1]) == 1.0
        @test MOI.get(model, MOI.VariablePrimalStart(), x[2]) === nothing
        Gurobi._update_if_necessary(model)
        @test Gurobi.get_dblattrelement(
            model.inner, "Start", Gurobi._info(model, x[2]).column
        ) == Gurobi.GRB_UNDEFINED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.ObjectiveValue()) == 0.0
        # We don't support ConstraintDualStart or ConstraintPrimalStart yet.
        # @test_broken MOIT.start_values_test(Gurobi.Optimizer(GUROBI_ENV), OPTIMIZER)
    end

    @testset "supports_constrainttest" begin
        # supports_constrainttest needs VectorOfVariables-in-Zeros,
        # MOIT.supports_constrainttest(Gurobi.Optimizer(GUROBI_ENV), Float64, Float32)
        # but supports_constrainttest is broken via bridges:
        MOI.empty!(OPTIMIZER)
        MOI.add_variable(OPTIMIZER)
        @test  MOI.supports_constraint(OPTIMIZER, MOI.SingleVariable, MOI.EqualTo{Float64})
        @test  MOI.supports_constraint(OPTIMIZER, MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})
        # This test is broken for some reason:
        @test_broken !MOI.supports_constraint(OPTIMIZER, MOI.ScalarAffineFunction{Int}, MOI.EqualTo{Float64})
        @test !MOI.supports_constraint(OPTIMIZER, MOI.ScalarAffineFunction{Int}, MOI.EqualTo{Int})
        @test !MOI.supports_constraint(OPTIMIZER, MOI.SingleVariable, MOI.EqualTo{Int})
        @test  MOI.supports_constraint(OPTIMIZER, MOI.VectorOfVariables, MOI.Zeros)
        @test !MOI.supports_constraint(OPTIMIZER, MOI.VectorOfVariables, MOI.EqualTo{Float64})
        @test !MOI.supports_constraint(OPTIMIZER, MOI.SingleVariable, MOI.Zeros)
        @test !MOI.supports_constraint(OPTIMIZER, MOI.VectorOfVariables, MOIT.UnknownVectorSet)
    end

    @testset "set_lower_bound_twice" begin
        MOIT.set_lower_bound_twice(OPTIMIZER, Float64)
    end

    @testset "set_upper_bound_twice" begin
        MOIT.set_upper_bound_twice(OPTIMIZER, Float64)
    end
end

@testset "LQOI Issue #38" begin
    # https://github.com/JuliaOpt/LinQuadOptInterface.jl/issues/38#issuecomment-407625187
    _getinner(opt::Gurobi.Optimizer) = opt.inner
    @inferred _getinner(Gurobi.Optimizer(GUROBI_ENV))
end

@testset "User limit handling (issue #140)" begin
    # Verify that we return the correct status codes when a mixed-integer
    # problem has been solved to a *feasible* but not necessarily optimal
    # solution. To do that, we will set up an intentionally dumbed-down
    # Gurobi Gurobi.Optimizer (with all heuristics and pre-solve turned off) and
    # ask it to solve a classic knapsack problem. Setting SolutionLimit=1
    # forces the solver to return after its first feasible MIP solution,
    # which tests the right part of the code without relying on potentially
    # flaky or system-dependent time limits.
    m = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0, SolutionLimit=1,
                         Heuristics=0.0, Presolve=0)
    N = 100
    x = MOI.add_variables(m, N)
    for xi in x
        MOI.add_constraint(m, MOI.SingleVariable(xi), MOI.ZeroOne())
        MOI.set(m, MOI.VariablePrimalStart(), xi, 0.0)
    end
    # Given a collection of items with individual weights and values,
    # maximize the total value carried subject to the constraint that
    # the total weight carried is less than 10.
    Random.seed!(1)
    item_weights = rand(N)
    item_values = rand(N)
    MOI.add_constraint(m,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0))
    MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(.-item_values, x), 0.0))
    MOI.optimize!(m)

    @test MOI.get(m, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    # We should have a primal feasible solution:
    @test MOI.get(m, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    # But we have no dual status:
    @test MOI.get(m, MOI.DualStatus()) == MOI.NO_SOLUTION
end

@testset "Constant objective (issue #111)" begin
    m = Gurobi.Optimizer(GUROBI_ENV)
    x = MOI.add_variable(m)
    MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 2.0))
    @test MOI.get(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()).constant == 2.0
    @test Gurobi.get_dblattr(m.inner, "ObjCon") == 2.0

    MOI.modify(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarConstantChange(3.0))
    @test MOI.get(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()).constant == 3.0
    @test Gurobi.get_dblattr(m.inner, "ObjCon") == 3.0
end

@testset "Env" begin
    @testset "User-provided" begin
        env = Gurobi.Env()
        model_1 = Gurobi.Optimizer(env)
        @test model_1.inner.env === env
        model_2 = Gurobi.Optimizer(env)
        @test model_2.inner.env === env
        # Check that finalizer doesn't touch env when manually provided.
        finalize(model_1.inner)
        @test Gurobi.is_valid(env)
    end
    @testset "Automatic" begin
        model_1 = Gurobi.Optimizer()
        model_2 = Gurobi.Optimizer()
        @test model_1.inner.env !== model_2.inner.env
        # Check that env is finalized with model when not supplied manually.
        finalize(model_1.inner)
        @test !Gurobi.is_valid(model_1.inner.env)
    end
    @testset "Env when emptied" begin
        @testset "User-provided" begin
            env = Gurobi.Env()
            model = Gurobi.Optimizer(env)
            @test model.inner.env === env
            @test Gurobi.is_valid(env)
            MOI.empty!(model)
            @test model.inner.env === env
            @test Gurobi.is_valid(env)
        end
        @testset "Automatic" begin
            model = Gurobi.Optimizer()
            env = model.inner.env
            MOI.empty!(model)
            @test model.inner.env !== env
            @test Gurobi.is_valid(model.inner.env)
        end
    end
end

@testset "Conflict refiner" begin
    @testset "Variable bounds (SingleVariable and LessThan/GreaterThan)" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(model)
        c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(2.0))
        c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
        @test_throws ErrorException MOI.get(model, Gurobi.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        Gurobi.compute_conflict(model)
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMAL
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c2) == true
    end

    @testset "Variable bounds (ScalarAffine)" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(model)
        c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(2.0))
        c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
        @test_throws ErrorException MOI.get(model, Gurobi.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        Gurobi.compute_conflict(model)
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMAL
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c2) == true
    end

    @testset "Variable bounds (Invali Interval)" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(model)
        c1 = MOI.add_constraint(
            model, MOI.SingleVariable(x), MOI.Interval(1.0, 0.0)
        )
        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
        @test_throws ErrorException MOI.get(model, Gurobi.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        Gurobi.compute_conflict(model)
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMAL
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c1) == true
    end

    @testset "Two conflicting constraints (GreaterThan, LessThan)" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
        b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
        cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
        c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
        cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
        c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
        @test_throws ErrorException MOI.get(model, Gurobi.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        Gurobi.compute_conflict(model)
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMAL
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), b1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), b2) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c2) == false
    end

    @testset "Two conflicting constraints (EqualTo)" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
        b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
        cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
        c1 = MOI.add_constraint(model, cf1, MOI.EqualTo(-1.0))
        cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
        c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
        @test_throws ErrorException MOI.get(model, Gurobi.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        Gurobi.compute_conflict(model)
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMAL
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), b1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), b2) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c2) == false
    end

    @testset "Variables outside conflict" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        z = MOI.add_variable(model)
        b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
        b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
        b3 = MOI.add_constraint(model, MOI.SingleVariable(z), MOI.GreaterThan(0.0))
        cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
        c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
        cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0, 1.0], [x, y, z]), 0.0)
        c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
        @test_throws ErrorException MOI.get(model, Gurobi.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        Gurobi.compute_conflict(model)
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMAL
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), b1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), b2) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), b3) == false
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c1) == true
        @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c2) == false
    end

    @testset "No conflict" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(model)
        c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(1.0))
        c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(2.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
        @test_throws ErrorException MOI.get(model, Gurobi.ConstraintConflictStatus(), c1)

        @test_throws Gurobi.GurobiError Gurobi.compute_conflict(model)
        # TODO(odow): bypass Gurobi's IIS checker when underlying model is
        # feasible.
        # @test MOI.get(model, Gurobi.ConflictStatus()) == MOI.INFEASIBLE
        # @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c1) == false
        # @test MOI.get(model, Gurobi.ConstraintConflictStatus(), c2) == false
    end
end

@testset "RawParameter" begin
    model = Gurobi.Optimizer(GUROBI_ENV)
    @test MOI.get(model, MOI.RawParameter("OutputFlag")) == 1
    MOI.set(model, MOI.RawParameter("OutputFlag"), 0)
    @test MOI.get(model, MOI.RawParameter("OutputFlag")) == 0
end

@testset "QCPDuals without needing to pass QCPDual=1" begin
    @testset "QCPDual=1" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0, QCPDual=1)
        MOI.Utilities.loadfromstring!(model, """
        variables: x, y, z
        minobjective: 1.0 * x + 1.0 * y + 1.0 * z
        c1: x + y == 2.0
        c2: x + y + z >= 0.0
        c3: 1.0 * x * x + -1.0 * y * y + -1.0 * z * z >= 0.0
        c4: x >= 0.0
        c5: y >= 0.0
        c6: z >= 0.0
        """)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
        c1 = MOI.get(model, MOI.ConstraintIndex, "c1")
        c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
        c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
        @test MOI.get(model, MOI.ConstraintDual(), c1) ≈ 1.0 atol=1e-6
        @test MOI.get(model, MOI.ConstraintDual(), c2) ≈ 0.0 atol=1e-6
        @test MOI.get(model, MOI.ConstraintDual(), c3) ≈ 0.0 atol=1e-6
    end
    @testset "QCPDual default" begin
        model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        MOI.Utilities.loadfromstring!(model, """
        variables: x, y, z
        minobjective: 1.0 * x + 1.0 * y + 1.0 * z
        c1: x + y == 2.0
        c2: x + y + z >= 0.0
        c3: 1.0 * x * x + -1.0 * y * y + -1.0 * z * z >= 0.0
        c4: x >= 0.0
        c5: y >= 0.0
        c6: z >= 0.0
        """)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION
        c1 = MOI.get(model, MOI.ConstraintIndex, "c1")
        c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
        c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
        @test_throws Gurobi.GurobiError MOI.get(model, MOI.ConstraintDual(), c1)
        @test_throws Gurobi.GurobiError MOI.get(model, MOI.ConstraintDual(), c2)
        @test_throws Gurobi.GurobiError MOI.get(model, MOI.ConstraintDual(), c3)
    end
end

@testset "Add and delete constraints" begin
    model = Gurobi.Optimizer(GUROBI_ENV)
    x = MOI.add_variables(model, 2)
    cs = MOI.add_constraints(
        model,
        [MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x[i])], 0.0) for i in 1:2],
        MOI.EqualTo.([0.0, 0.0])
    )
    @test MOI.get(model, MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}
    }()) == 2
    MOI.delete(model, cs)
    @test iszero(MOI.get(model, MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}
    }()))
end

@testset "Extra name tests" begin
    model = Gurobi.Optimizer(GUROBI_ENV)
    @testset "Variables" begin
        MOI.empty!(model)
        x = MOI.add_variables(model, 3)
        MOI.set(model, MOI.VariableName(), x[1], "x1")
        @test MOI.get(model, MOI.VariableIndex, "x1") == x[1]
        MOI.set(model, MOI.VariableName(), x[1], "x2")
        @test MOI.get(model, MOI.VariableIndex, "x1") === nothing
        @test MOI.get(model, MOI.VariableIndex, "x2") == x[1]
        MOI.set(model, MOI.VariableName(), x[2], "x1")
        @test MOI.get(model, MOI.VariableIndex, "x1") == x[2]
        MOI.set(model, MOI.VariableName(), x[3], "xα")
        @test MOI.get(model, MOI.VariableIndex, "xα") == x[3]
        MOI.set(model, MOI.VariableName(), x[1], "x1")
        @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x1")
    end

    @testset "Variable bounds" begin
        MOI.empty!(model)
        x = MOI.add_variable(model)
        c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
        c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))
        MOI.set(model, MOI.ConstraintName(), c1, "c1")
        @test MOI.get(model, MOI.ConstraintIndex, "c1") == c1
        MOI.set(model, MOI.ConstraintName(), c1, "c2")
        @test MOI.get(model, MOI.ConstraintIndex, "c1") === nothing
        @test MOI.get(model, MOI.ConstraintIndex, "c2") == c1
        MOI.set(model, MOI.ConstraintName(), c2, "c1")
        @test MOI.get(model, MOI.ConstraintIndex, "c1") == c2
        MOI.set(model, MOI.ConstraintName(), c1, "c1")
        @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "c1")
    end

    @testset "Affine constraints" begin
        MOI.empty!(model)
        x = MOI.add_variable(model)
        f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
        c1 = MOI.add_constraint(model, f, MOI.GreaterThan(0.0))
        c2 = MOI.add_constraint(model, f, MOI.LessThan(1.0))
        MOI.set(model, MOI.ConstraintName(), c1, "c1")
        @test MOI.get(model, MOI.ConstraintIndex, "c1") == c1
        MOI.set(model, MOI.ConstraintName(), c1, "c2")
        @test MOI.get(model, MOI.ConstraintIndex, "c1") === nothing
        @test MOI.get(model, MOI.ConstraintIndex, "c2") == c1
        MOI.set(model, MOI.ConstraintName(), c2, "c1")
        @test MOI.get(model, MOI.ConstraintIndex, "c1") == c2
        MOI.set(model, MOI.ConstraintName(), c1, "c1")
        @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "c1")
    end
end

@testset "ConstraintAttribute" begin
    model = Gurobi.Optimizer(GUROBI_ENV)
    MOI.Utilities.loadfromstring!(model, """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""")
    c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
    c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
    # Linear constraints are supported - one test for each different type.
    # Integer attribute
    MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c2, 2)
    @test MOI.get(model, Gurobi.ConstraintAttribute("Lazy"), c2) == 2
    # Real attribute
    MOI.set(model, Gurobi.ConstraintAttribute("RHS"), c2, 2.0)
    @test MOI.get(model, Gurobi.ConstraintAttribute("RHS"), c2) == 2.0
    # Char attribute
    MOI.set(model, Gurobi.ConstraintAttribute("Sense"), c2, '<')
    @test MOI.get(model, Gurobi.ConstraintAttribute("Sense"), c2) == '<'
    # String attribute
    MOI.set(model, Gurobi.ConstraintAttribute("ConstrName"), c2, "c4")
    @test MOI.get(model, Gurobi.ConstraintAttribute("ConstrName"), c2) == "c4"
    # Things that should fail follow.
    # Non-linear constraints are not supported.
    @test_throws(
        MOI.SetAttributeNotAllowed(Gurobi.ConstraintAttribute("Lazy")),
        MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c3, 1)
    )
    # Getting/setting a non-existing attribute.
    attr = Gurobi.ConstraintAttribute("Non-existing")
    @test_throws MOI.UnsupportedAttribute(attr) MOI.set(model, attr, c2, 1)
    @test_throws MOI.UnsupportedAttribute(attr) MOI.get(model, attr, c2)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError("Attribute Lazy is Integer but Float64 provided."),
        MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c2, 1.0)
    )
end

@testset "VariableAttribute" begin
    model = Gurobi.Optimizer(GUROBI_ENV)
    MOI.Utilities.loadfromstring!(model, """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""")
    x = MOI.get(model, MOI.VariableIndex, "x")
    # Setting attributes of each type
    # Integer attribute
    MOI.set(model, Gurobi.VariableAttribute("VarHintPri"), x, 2)
    @test MOI.get(model, Gurobi.VariableAttribute("VarHintPri"), x) == 2
    # Real Attribute
    MOI.set(model, Gurobi.VariableAttribute("LB"), x, 2.0)
    @test MOI.get(model, Gurobi.VariableAttribute("LB"), x) == 2.0
    # Char Attribute
    MOI.set(model, Gurobi.VariableAttribute("VType"), x, 'B')
    @test MOI.get(model, Gurobi.VariableAttribute("VType"), x) == 'B'
    # String Attribute
    MOI.set(model, Gurobi.VariableAttribute("VarName"), x, "my_var")
    @test MOI.get(model, Gurobi.VariableAttribute("VarName"), x) == "my_var"
    # Things that should fail follow.
    # Getting/setting a non-existing attribute.
    attr = Gurobi.VariableAttribute("Non-existing")
    @test_throws MOI.UnsupportedAttribute(attr) MOI.set(model, attr, x, 1)
    @test_throws MOI.UnsupportedAttribute(attr) MOI.get(model, attr, x)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError("Attribute BranchPriority is Integer but Float64 provided."),
        MOI.set(model, Gurobi.VariableAttribute("BranchPriority"), x, 1.0)
    )
end

@testset "ModelAttribute" begin
    model = Gurobi.Optimizer(GUROBI_ENV)
    MOI.Utilities.loadfromstring!(model, """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""")
    # Setting attributes of each type
    # Integer attribute
    MOI.set(model, Gurobi.ModelAttribute("ModelSense"), -1)
    @test MOI.get(model, Gurobi.ModelAttribute("ModelSense")) == -1
    # Real Attribute
    MOI.set(model, Gurobi.ModelAttribute("ObjCon"), 3.0)
    @test MOI.get(model, Gurobi.ModelAttribute("ObjCon")) == 3.0
    # String Attribute
    MOI.set(model, Gurobi.ModelAttribute("ModelName"), "My model")
    @test MOI.get(model, Gurobi.ModelAttribute("ModelName")) == "My model"
    # Things that should fail follow.
    # Getting/setting a non-existing attribute.
    attr = Gurobi.ModelAttribute("Non-existing")
    @test_throws MOI.UnsupportedAttribute(attr) MOI.set(model, attr, 1)
    @test_throws MOI.UnsupportedAttribute(attr) MOI.get(model, attr)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError("Attribute NumStart is Integer but Float64 provided."),
        MOI.set(model, Gurobi.ModelAttribute("NumStart"), 4.0)
    )
end

@testset "SOC hide lower bound constraint" begin
    model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
    @testset "No initial bound" begin
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        MOI.add_constraints(
            model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
        )
        c_soc = MOI.add_constraint(
            model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
        )
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_soc)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
    end
    @testset "non-negative initial bound" begin
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(1.0))
        MOI.add_constraints(
            model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
        )
        c_soc = MOI.add_constraint(
            model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
        )
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_soc)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 1.0
    end
    @testset "negative initial bound" begin
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-1.0))
        MOI.add_constraints(
            model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
        )
        c_soc = MOI.add_constraint(
            model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
        )
        MOI.optimize!(model)
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_soc)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == -1.0
    end
    @testset "non-negative post bound" begin
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        MOI.add_constraints(
            model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
        )
        MOI.add_constraint(
            model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
        )
        c_lb = MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(6.0))
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 6.0
        MOI.delete(model, c_lb)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    end
    @testset "negative post bound" begin
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        MOI.add_constraints(
            model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
        )
        c_soc = MOI.add_constraint(
            model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
        )
        MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-6.0))
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_soc)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == -6.0
    end
    @testset "negative post bound II" begin
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        MOI.add_constraints(
            model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
        )
        c_soc = MOI.add_constraint(
            model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
        )
        c_lb = MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-6.0))
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_lb)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_soc)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
    end
    @testset "negative post bound III" begin
        # This test was added because add_constraint and add_constraints had
        # different implementations and add_constraints failed where
        # add_constraint succeeded.
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        c_soc = MOI.add_constraint(
            model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
        )
        c_lbs = MOI.add_constraints(
            model, MOI.SingleVariable.([t; x]), MOI.GreaterThan.([-6.0, 3.0, 4.0])
        )
        c_lb = first(c_lbs)
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_lb)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
        MOI.delete(model, c_soc)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
    end
    @testset "2 SOC's" begin
        MOI.empty!(model)
        t = MOI.add_variable(model)
        x = MOI.add_variables(model, 2)
        MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-1.0))
        MOI.add_constraints(
            model, MOI.SingleVariable.(x), MOI.GreaterThan.([4.0, 3.0])
        )
        c_soc_1 = MOI.add_constraint(
            model, MOI.VectorOfVariables([t, x[1]]), MOI.SecondOrderCone(2)
        )
        c_soc_2 = MOI.add_constraint(
            model, MOI.VectorOfVariables([t, x[2]]), MOI.SecondOrderCone(2)
        )
        MOI.optimize!(model)
        MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 4.0
        MOI.delete(model, c_soc_1)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == 3.0
        MOI.delete(model, c_soc_2)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), t) == -1.0
    end
end

@testset "Duplicate names" begin
    @testset "Variables" begin
        model = Gurobi.Optimizer(GUROBI_ENV)
        (x, y, z) = MOI.add_variables(model, 3)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "x")
        MOI.set(model, MOI.VariableName(), z, "z")
        @test MOI.get(model, MOI.VariableIndex, "z") == z
        @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        @test MOI.get(model, MOI.VariableIndex, "x") == x
        @test MOI.get(model, MOI.VariableIndex, "y") == y
        MOI.set(model, MOI.VariableName(), z, "x")
        @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x")
        MOI.delete(model, x)
        @test MOI.get(model, MOI.VariableIndex, "x") == z
    end
    @testset "SingleVariable" begin
        model = Gurobi.Optimizer(GUROBI_ENV)
        x = MOI.add_variables(model, 3)
        c = MOI.add_constraints(model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0))
        MOI.set(model, MOI.ConstraintName(), c[1], "x")
        MOI.set(model, MOI.ConstraintName(), c[2], "x")
        MOI.set(model, MOI.ConstraintName(), c[3], "z")
        @test MOI.get(model, MOI.ConstraintIndex, "z") == c[3]
        @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
        MOI.set(model, MOI.ConstraintName(), c[2], "y")
        @test MOI.get(model, MOI.ConstraintIndex, "x") == c[1]
        @test MOI.get(model, MOI.ConstraintIndex, "y") == c[2]
        MOI.set(model, MOI.ConstraintName(), c[3], "x")
        @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
        MOI.delete(model, c[1])
        @test MOI.get(model, MOI.ConstraintIndex, "x") == c[3]
        MOI.set(model, MOI.ConstraintName(), c[2], "x")
        @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
        MOI.delete(model, x[3])
        @test MOI.get(model, MOI.ConstraintIndex, "x") == c[2]
    end
    @testset "ScalarAffineFunction" begin
        model = Gurobi.Optimizer(GUROBI_ENV)
        x = MOI.add_variables(model, 3)
        fs = [
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, xi)], 0.0)
            for xi in x
        ]
        c = MOI.add_constraints(model, fs, MOI.GreaterThan(0.0))
        MOI.set(model, MOI.ConstraintName(), c[1], "x")
        MOI.set(model, MOI.ConstraintName(), c[2], "x")
        MOI.set(model, MOI.ConstraintName(), c[3], "z")
        @test MOI.get(model, MOI.ConstraintIndex, "z") == c[3]
        @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
        MOI.set(model, MOI.ConstraintName(), c[2], "y")
        @test MOI.get(model, MOI.ConstraintIndex, "x") == c[1]
        @test MOI.get(model, MOI.ConstraintIndex, "y") == c[2]
        MOI.set(model, MOI.ConstraintName(), c[3], "x")
        @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
        MOI.delete(model, c[1])
        @test MOI.get(model, MOI.ConstraintIndex, "x") == c[3]
    end
end

@testset "Duals with equal bounds #250" begin
    model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
    x = MOI.add_variable(model)
    xl = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    xu = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ConstraintDual(), xl) == 1.0
    @test MOI.get(model, MOI.ConstraintDual(), xu) == 0.0
end

@testset "Objective functions" begin
    model = Gurobi.Optimizer(GUROBI_ENV)
    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    @test MOI.get(model, MOI.ListOfModelAttributesSet()) == Any[MOI.ObjectiveSense()]
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))
    @test MOI.get(model, MOI.ListOfModelAttributesSet()) ==
        Any[MOI.ObjectiveSense(), MOI.ObjectiveFunction{MOI.SingleVariable}()]
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(model, MOI.ListOfModelAttributesSet()) == Any[MOI.ObjectiveSense()]
end

@testset "FEASIBILITY_SENSE zeros objective" begin
    model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) == 1.0
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) == 0.0
end

@testset "Attributes" begin
    model = Gurobi.Optimizer(GUROBI_ENV, OutputFlag = 0)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.NodeCount()) == 0
    @test MOI.get(model, MOI.BarrierIterations()) == 0
    @test MOI.get(model, MOI.SimplexIterations()) == 0
end


text:
macro err_check(expr)
    quote
        ret = $expr
        ret == 0 || mbed_err(ret)
        ret
    end
end

immutable MbedException <: Exception
    ret::Cint
end

function show(io::IO, err::MbedException)
    print(io, "MbedTLS error code $(err.ret): $(strerror(err.ret))")
end

mbed_err(ret) = throw(MbedException(ret))

function strerror(ret, bufsize=1000)
    buf = Vector{UInt8}(bufsize)
    ccall((:mbedtls_strerror, MBED_CRYPTO), Cint,
        (Cint, Ptr{Void}, Csize_t),
        ret, buf, bufsize)
    String(pointer(buf))
end

ext: jl


text:
using NFFT, DataFrames, LinearAlgebra, LaTeXStrings, DelimitedFiles
using BenchmarkTools
using Plots; pgfplotsx()
using Plots.Measures

include("../../Wrappers/NFFT3.jl")
include("../../Wrappers/FINUFFT.jl")

const packagesCtor = [NFFTPlan, NFFTPlan, NFFT3Plan, FINUFFTPlan] #, NFFTPlan
const packagesStr = [ "NFFT.jl/TENSOR", "NFFT.jl/POLY", "NFFT3/TENSOR", "FINUFFT"] #"NFFT.jl/LINEAR"
const precomp = [NFFT.TENSOR, NFFT.POLYNOMIAL, NFFT.TENSOR, NFFT.LINEAR] #NFFT.LINEAR
const blocking = [true, true, true, true, true]

const benchmarkTime = [20, 80]

NFFT.FFTW.set_num_threads(Threads.nthreads())
ccall(("omp_set_num_threads",NFFT3.lib_path_nfft),Nothing,(Int64,),convert(Int64,Threads.nthreads()))
@info ccall(("nfft_get_num_threads",NFFT3.lib_path_nfft),Int64,())
NFFT._use_threads[] = (Threads.nthreads() > 1)

const σs = [2.0]
const ms = 3:8
#const NBase = [65536, 256, 32] #[4*4096, 128, 32]
#const NBase = [4*4096, 128, 32]
const NBase = [512*512, 512, 64]
const Ds = 1:3
const fftflags = NFFT.FFTW.MEASURE

function nfft_accuracy_comparison(Ds=1:3)
  println("\n\n ##### nfft_performance vs accuracy ##### \n\n")

  df = DataFrame(Package=String[], D=Int[], M=Int[], N=Int[], m = Int[], σ=Float64[],
                   ErrorTrafo=Float64[], ErrorAdjoint=Float64[],
                   TimePre=Float64[], TimeTrafo=Float64[], TimeAdjoint=Float64[],  )

  for D in Ds
    @info "### Dimension D=$D ###"
    N = ntuple(d->NBase[D], D)
    M = prod(N)

    x = rand(D,M) .- 0.5
    fHat = randn(ComplexF64, M)
    fApprox = randn(ComplexF64, N)
    gHatApprox = randn(ComplexF64, M)

    # ground truth (numerical)
    pNDFT = NDFTPlan(x, N)
    f = adjoint(pNDFT) * fHat
    gHat = pNDFT * f

    for σ in σs
      for m in ms
        @info "m=$m D=$D σ=$σ "

        for pl = 1:length(packagesStr)
          planner = packagesCtor[pl]
          p = planner(x, N; m, σ, precompute=precomp[pl], blocking=blocking[pl], fftflags=fftflags)
          b = @benchmark $planner($x, $N; m=$m, σ=$σ, precompute=$(precomp[pl]), blocking=$(blocking[pl]), fftflags=$(fftflags))
          tpre = minimum(b).time / 1e9

          @info "Adjoint accuracy: $(packagesStr[pl])"
          mul!(fApprox, adjoint(p), fHat)
          eadjoint = norm(f[:] - fApprox[:]) / norm(f[:])

          @info "Adjoint benchmark: $(packagesStr[pl])"
          BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[1]
          b = @benchmark mul!($fApprox, $(adjoint(p)), $fHat)
          tadjoint = minimum(b).time / 1e9

          @info "Trafo accuracy: $(packagesStr[pl])"
          mul!(gHatApprox, p, f)
          etrafo = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])

          @info "Trafo benchmark: $(packagesStr[pl])"
          BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[2]
          b = @benchmark mul!($gHatApprox, $p, $f)
          ttrafo = minimum(b).time / 1e9

          push!(df, (packagesStr[pl], D, M, N[D], m, σ, etrafo, eadjoint, tpre, ttrafo, tadjoint))

        end
      end
    end
  end
  return df
end

function plot_accuracy(df, packagesStr, packagesStrShort, filename)

  Plots.scalefontsizes()
  Plots.scalefontsizes(1.5)

  #colors = [RGB(0.0,0.29,0.57), RGB(0.3,0.5,0.7), RGB(0.95,0.59,0.22), RGB(1.0,0.87,0.0)]
  #colors = [RGB(0.0,0.29,0.57), RGB(0.3,0.5,0.7),  RGB(0.7,0.13,0.16), RGB(0.72,0.84,0.48)]
  #colors = [RGB(0.0,0.29,0.57), RGB(0.3,0.5,0.7), RGB(1.0,0.87,0.0), RGB(0.95,0.59,0.22)]
  colors = [RGB(0.0,0.29,0.57), RGB(0.3,0.5,0.7), RGB(0.94,0.53,0.12), RGB(0.99,0.75,0.05)]
  ls = [:solid, :solid, :solid, :solid]
  shape = [:xcross, :circle, :xcross, :cross]

  xlims = [(4e-13,1e-5), (4e-15,1e-4),(4e-15,1e-4)]

  pl = Matrix{Any}(undef, 3, length(Ds))
  for (i,D) in enumerate(Ds)
    titleTrafo = L"\textrm{NFFT}, \textrm{%$(D)D}"
    titleAdjoint = L"\textrm{NFFT}^H, \textrm{%$(D)D}"
    titlePre = L"\textrm{Precompute}, \textrm{%$(D)D}"
    xlabel = (i==length(Ds)) ? "Relative Error" : ""

    df1_ = df[df.σ.==2.0 .&& df.D.==D,:]
    maxTimeTrafo = maximum( maximum(df1_[df1_.Package.==pStr,:TimeTrafo]) for pStr in packagesStr)
    maxTimeAdjoint = maximum( maximum(df1_[df1_.Package.==pStr,:TimeAdjoint]) for pStr in packagesStr)
    maxTimePre = maximum( maximum(df1_[df1_.Package.==pStr,:TimePre]) for pStr in packagesStr)

    p1 = plot(df1_[df1_.Package.==packagesStr[1],:ErrorTrafo],
              df1_[df1_.Package.==packagesStr[1],:TimeTrafo], ylims=(0.0,maxTimeTrafo),
              label = packagesStrShort[1],
              xscale = :log10, legend = (i==length(Ds)) ? (0.0, -0.5) : nothing,
              lw=2, xlabel = xlabel, ylabel="Runtime / s",
              title=titleTrafo, shape=shape[1], ls=ls[1],
              c=colors[1], msc=colors[1], mc=colors[1], ms=4, msw=2,
              xlims=xlims[i])

    for p=2:length(packagesStr)
      plot!(p1, df1_[df1_.Package.==packagesStr[p],:ErrorTrafo],
            df1_[df1_.Package.==packagesStr[p],:TimeTrafo],
              xscale = :log10, lw=2, shape=shape[p], ls=ls[p],
              label =  packagesStrShort[p] ,
              c=colors[p], msc=colors[p], mc=colors[p], ms=4, msw=2)
    end

    p2 = plot(df1_[df1_.Package.==packagesStr[1],:ErrorAdjoint],
              df1_[df1_.Package.==packagesStr[1],:TimeAdjoint], ylims=(0.0,maxTimeAdjoint),
              xscale = :log10,  lw=2, xlabel = xlabel, #ylabel="Runtime / s", #label=packagesStr[1],
              legend = nothing, title=titleAdjoint, shape=shape[1], ls=ls[1],
              c=colors[1], msc=colors[1], mc=colors[1], ms=4, msw=2,
              xlims=xlims[i])

    for p=2:length(packagesStr)
      plot!(p2, df1_[df1_.Package.==packagesStr[p],:ErrorAdjoint],
            df1_[df1_.Package.==packagesStr[p],:TimeAdjoint],
              xscale = :log10,  lw=2, shape=shape[p], ls=ls[p], #label=packagesStr[p],
              c=colors[p], msc=colors[p], mc=colors[p], ms=4, msw=2)
    end

    p3 = plot(df1_[df1_.Package.==packagesStr[1],:ErrorAdjoint],
              df1_[df1_.Package.==packagesStr[1],:TimePre], ylims=(0.0,maxTimePre),
              label = (i==1) ? packagesStrShort[1] : "",
              xscale = :log10,  lw=2, xlabel = xlabel, #ylabel="Runtime / s", #label=packagesStr[1],
              title=titlePre, shape=shape[1], ls=ls[1],
              c=colors[1], msc=colors[1], mc=colors[1], ms=4, msw=2,
              legend = nothing,
              xlims=xlims[i] )

    for p=2:length(packagesStr)
      plot!(p3, df1_[df1_.Package.==packagesStr[p],:ErrorAdjoint],
            df1_[df1_.Package.==packagesStr[p],:TimePre],
             label = (i==1) ? packagesStrShort[p] : "",
              xscale = :log10,  lw=2, shape=shape[p], ls=ls[p], #label=packagesStr[p],
              c=colors[p], msc=colors[p], mc=colors[p], ms=4, msw=2)
    end
    pl[1,i] = p1; pl[2,i] = p2; pl[3,i] = p3;
  end

  p = plot(vec(pl)..., layout=(length(Ds),3), size=(900,600), dpi=200, margin = 1mm)

  mkpath("./img/")
  savefig(p, filename)
  return p
end

#df = nfft_accuracy_comparison(Ds)
#writedlm("data/performanceVsAccuracy.csv", Iterators.flatten(([names(df)], eachrow(df))), ',')

data, header = readdlm("data/performanceVsAccuracy.csv", ',', header=true);
df = DataFrame(data, vec(header))

plot_accuracy(df, [ "NFFT.jl/POLY", "NFFT.jl/TENSOR", "NFFT3/TENSOR", "FINUFFT"],
                  [ "NFFT.jl/POLY", "NFFT.jl/TENSOR", "NFFT3", "FINUFFT"], "./img/performanceVsAccuracy.pdf")

#plot_accuracy(df, [ "NFFT.jl/POLY", "NFFT.jl/TENSOR" ], #"NFFT.jl/LINEAR"  , "LINEAR"
#                  [ "POLYNOMIAL", "TENSOR"], "./img/performanceVsAccuracyPrecomp.pdf")


text:
using BNPRegressionGGA2021
using CSV
using DataFrames
using Random
const BNP = BNPRegressionGGA2021

# Simulate the data
begin
    include("dgp-normal-normal2.jl");
    Random.seed!(1);
    N0, N1 = 500, 50;
    dy, y0, x0, X0, y1, x1, X1 = simulate_sample(N0, N1);
    y0 = (y0 .- mean(y0)) ./ std(y0)
end;

# Fit the model
begin
    update_γ = false
    mapping = [[1], collect(2:size(X0, 2))];
    m = BNP.DGSBPNormal(; y0, X0, y1, X1, mapping, update_γ);
    chainf, chainβ, chaing = BNP.sample!(m; mcmcsize = 25000, burnin = 5000);
end;

# Save the results as a CSV file
begin
    df = DataFrame(
        y = y1,
        x = x1,
        fh = mean(chainf),
        # f0 = pdf.(dy.(x1), y1)
    )
    filename = "data/dgp-normal-normal-example2.csv"
    CSV.write(filename, df)
end


text:
"""

"""
function potentialflow(ctx, A, V, W)
    n = getproperty(ctx, :n)

    Φ = scalar(:Φ, n)

    U = gradient(neu, ctx, A, V, 0, Φ, 0)
    ΔΦ = divergence(neu, ctx, A, U, W)

    sym = linearize(star, ctx, ΔΦ, Φ)

    # linearity of gradient and divergence
    num = Dict{Int,Vector{Float64}}()
    for (key, val) in sym
        num[key] = getproperty.(val, Ref(:val))
    end

    Δ = spdiagm(num...)

    for i in axes(Δ, 1)
        Δ[i, i] = abs(Δ[i, i]) ≤ eps(Δ[i, i]) ? 1. : Δ[i, i]
    end

    U = gradient(neu, ctx, A, V, 0, zeros(prod(n)), 0)
    b = divergence(neu, ctx, A, U, W)

    phi = -Δ \ b

    phi, gradient(neu, ctx, A, V, 0, phi, 0)
end


text:
using BHAPtfem
using Test

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  @test materials[:steel].sg == 0.283

end # cd


text:
println("Testing floating normal progress bar (offset 4)")
function testfunc1(n, dt, tsleep, desc, barlen, offset)
    p = ProgressMeter.Progress(n, dt, desc, barlen, offset=offset)
    for i = 1:n
        sleep(tsleep)
        ProgressMeter.next!(p)
    end
    print("\n" ^ 5)
end
testfunc1(50, 0.2, 0.2, "progress  ", 70, 4)

println("Testing floating normal progress bars with values and keep (2 levels)")
function testfunc2(n, dt, tsleep, desc, barlen)
    p1 = ProgressMeter.Progress(n, dt, desc, barlen, offset=0)
    p2 = ProgressMeter.Progress(n, dt, desc, barlen, offset=5)
    for i = 1:n
        sleep(tsleep)
        ProgressMeter.next!(p1; showvalues = [(:i, i),
            (:constant, "foo"), (:isq, i^2), (:large, 2^i)], keep=false)
        ProgressMeter.next!(p2; showvalues = [(:i, i),
            (:constant, "foo"), (:isq, i^2), (:large, 2^i)])
    end
    print("\n" ^ 10)
end
testfunc2(50, 0.2, 0.2, "progress  ", 70)

println("Testing floating normal progress bars with changing offset")
function testfunc3(n, dt, tsleep, desc, barlen)
    p1 = ProgressMeter.Progress(n, dt, desc, barlen, offset=0)
    p2 = ProgressMeter.Progress(n, dt, desc, barlen, offset=1)
    for i = 1:n
        sleep(tsleep)
        ProgressMeter.next!(p1; showvalues = [(:i, i) for _ in 1:i], keep=false)
        ProgressMeter.next!(p2; showvalues = [(:i, i),
            (:constant, "foo"), (:isq, i^2), (:large, 2^i)], offset = (p1.offset + p1.numprintedvalues))
    end
    print("\n" ^ (10 + 5))
end
testfunc3(10, 0.2, 0.5, "progress  ", 70)

println("Testing floating thresh progress bar (offset 2)")
function testfunc4(thresh, dt, tsleep, desc, offset)
    prog = ProgressMeter.ProgressThresh(thresh, dt, desc, offset=offset)
    for val in 10 .^ range(2, stop=-6, length=20)
        ProgressMeter.update!(prog, val)
        sleep(tsleep)
    end
    print("\n" ^ 3)
end
testfunc4(1e-5, 0.2, 0.2, "Minimizing: ", 2)

println("Testing floating in @showprogress macro (3 levels)")
function testfunc5(n, tsleep)
    ProgressMeter.@showprogress "Level 0 " for i in 1:n
        ProgressMeter.@showprogress " Level 1 " 1 for i2 in 1:n
            ProgressMeter.@showprogress "  Level 2 " 2 for i3 in 1:n
                sleep(tsleep)
            end
        end
    end
    print("\n" ^ 2)
end
testfunc5(5, 0.1)


text:
# identity.jl

import LinearAlgebra: UniformScaling

"""
    *(I, X) = X
    *(J, X) = J.λ * X

Extends the effect of `I::UniformScaling` and scaled versions thereof
to also apply to `X::AbstractArray` instead of just to `Vector` and `Matrix` types.
"""
Base.:(*)(J::UniformScaling, X::AbstractArray) = J.λ == 1 ? X : J.λ * X


text:
# Default Color Schemes
# catergory 20b and 20c used as reference
const color_schemes=Dict{Symbol,Any}(
    :blues => colorscheme2array(ColorSchemes.ColorScheme(range(Colors.colorant"#3182BD", Colors.colorant"#C6DBEF", length=5))),#["#$(Colors.hex(c))" for c in ColorSchemes.ColorScheme(range(Colors.colorant"#3182BD", Colors.colorant"#C6DBEF", length=5))],
    :greens => colorscheme2array(ColorSchemes.ColorScheme(range(Colors.colorant"#31A354", Colors.colorant"#C7E9C0", length=5))),
    :oranges => colorscheme2array(ColorSchemes.ColorScheme(range(Colors.colorant"#E6550D", Colors.colorant"#FDD0A2", length=5))),
    :reds => colorscheme2array(ColorSchemes.ColorScheme(range(Colors.colorant"#843C39", Colors.colorant"#E7969C", length=5))),
    :purples => colorscheme2array(ColorSchemes.ColorScheme(range(Colors.colorant"#756BB1", Colors.colorant"#DADAEB", length=5))),

)

# Default plot attributes
default_plot_attributes = Dict{Symbol, Any}(
  :gen_color => color_schemes[:oranges],
  :bus_color => color_schemes[:greens],
  :branch_color => color_schemes[:blues],
  :connector_color => [:gray],
  :dcline_color => color_schemes[:purples],
  :storage_color => color_schemes[:oranges],
  :gen_size => 5e2,
  :bus_size => 5e2,
  :branch_size => 5,
  :connector_size => 3,
  :dcline_size => 5,
  :storage_size => 5e2,
  :width => 500,
  :height => 500,
  :gen_data => "ComponentType",
  :bus_data => "ComponentType",
  :branch_data => "ComponentType",
  :dcline_data => "ComponentType",
  :storage_data => "ComponentType",
  :gen_data_type => "nominal",
  :bus_data_type => "nominal",
  :branch_data_type => "nominal",
  :dcline_data_type => "nominal",
  :storage_data_type => "nominal",
);

# Returns a deepcopy of the default_plot_attributes, used to initialize plot_attributes in utils.jl
function copy_default_attributes()
  return deepcopy(default_plot_attributes);
end

# Dictionary of aliases of properties
# Aliases are replaced by their key in plot_attributes
# for example, nodecolor/markercolor/marker_color, edgecolor/edge_color/ec can be used to
# set nodecolor and edgecolor respectively
const attribute_aliases = Dict(
  :gen_color => [:gencolor, :generatorcolor, :generator_color, :node_color],
  :bus_color => [:buscolor, :substationcolor, :substation_color, :node_color],
  :branch_color => [:branchcolor, :transmissionlinecolor, :transmissionline_color,:transmission_line_color, :edge_color],
  :connector_color => [:connectorcolor, :edge_color],
  :storage_color => [:storagecolor, :batterycolor, :battery_color, :node_color],
  :dcline_color => [:dclinecolor, :dc_line_color, :edge_color],
  :bus_size => [:bussize, :node_size],
  :gen_size => [:gensize, :node_size],
  :storage_size => [:storagesize, :node_size],
  :branch_size => [:branchsize, :line_size, :edge_size],
  :dcline_size => [:dclinesize, :line_size, :edge_size],
  :connector_size => [:connectorsize, :edge_size],
)

const _color_attributes = [ # color (String or Symbol) type parameters
  :gen_color,
  :bus_color,
  :branch_color,
  :connector_color,
  :storage_color,
  :dcline_color
]
const _numeric_attributes = [ # numeric parameters
  :gen_size,
  :bus_size,
  :branch_size,
  :connector_size,
  :dcline_size,
  :storage_size,
  :width,
  :height
]
const _label_attributes = [ # label (String or Symbol) type parameters
  :gen_data,
  :bus_data,
  :branch_data,
  :dcline_data,
  :storage_data,
  :gen_data_type,
  :bus_data_type,
  :branch_data_type,
  :dcline_data_type,
  :storage_data_type
]


text:
#=
WARNING: the other schedulers are mostly experimental and likely to disappear in
the future.
=#

"""
    struct RunnableChannel <: AbstractChannel{DataFlowTask}

Used to store tasks which have been tagged as dependency-free, and thus can be
executed. The underlying data is stored using a priority queue, with elements
with a high priority being popped first.

Calling `take` on an empty `RunnableChannel` will block.
"""
struct RunnableChannel{T}
    data::PriorityQueue{T,Float64}
    cond_take::Threads.Condition
    function RunnableChannel{T}() where {T}
        lock = Threads.ReentrantLock()
        cond_take = Threads.Condition(lock)
        # larger is more urgent in priority
        data = PriorityQueue{T,Float64}(Base.Order.ReverseOrdering())
        new(data,cond_take)
    end
end
Base.lock(c::RunnableChannel)   = lock(c.cond_take)
Base.unlock(c::RunnableChannel) = unlock(c.cond_take)

function Base.take!(c::RunnableChannel)
    lock(c)
    try
        while isempty(c.data)
            wait(c.cond_take)
        end
        @trace "length_runnable $(length(c.data)) $(time_ns())"
        v = dequeue!(c.data)
        @trace "length_runnable $(length(c.data)) $(time_ns())"
        return v
    finally
        unlock(c)
    end
end

function Base.put!(c::RunnableChannel{T},t::T,p=1) where {T}
    lock(c)
    try
        @trace "length_runnable $(length(c.data)) $(time_ns())"
        push!(c.data, t=>p)
        @trace "length_runnable $(length(c.data)) $(time_ns())"
        notify(c.cond_take)
    finally
        unlock(c)
    end
    return t
end

"""
    StaticScheduler{T} <: TaskGraphScheduler{T}

Like the [`JuliaScheduler`](@ref), but requires an explicit call to [`execute_dag(ex)`](@ref) to
start running the nodes in its `dag` (and removing them as they are completed).

Using a `StaticScheduler` is useful if you wish examine the
underling `TaskGraph` before it is executed.
"""
struct StaticScheduler{T} <: TaskGraphScheduler
    dag::DAG{T}
    finished::FinishedChannel{T}
    dag_worker::Task
    function StaticScheduler{T}() where {T}
        dag            = DAG{T}() # unbuffered dag
        finished       = FinishedChannel{T}()
        worker         = consume_finished(dag,finished)
        sch            = new(dag,finished,worker)
        return sch
    end
end
StaticScheduler() = StaticScheduler{DataFlowTask}()

# interface methods for StaticScheduler
function spawn(tj,::StaticScheduler)
    tj.task.sticky = false
    return tj
end

function Base.schedule(tj,::StaticScheduler)
    # delay the scheduling of tj.task until execute_dag is called
    return tj
end

"""
    execute_dag(sch::StaticScheduler)

Execute all the nodes in the task graph, removing them from the `dag` as they
are completed. This function waits for the dag to be emptied before returning.
"""
function execute_dag(sch::StaticScheduler)
    if isempty(sch.dag)
        @warn "DAG is empty"
        return sch
    else
        for (t,_) in sch.dag
            schedule(t.task)
        end
        # wait on dag to be emptied. This means all nodes have been ran and finished.
        wait(sch)
        return sch
    end
end

"""
    struct PriorityScheduler{T} <: TaskGraphScheduler{T}

Execute a `DAG` by spawning workers that take elements from the `runnable`
channel, execute them, and put them into a `finished` channel to be processed by
a `dag_worker`.
"""
struct PriorityScheduler{T} <: TaskGraphScheduler
    dag::DAG{T}
    runnable::RunnableChannel{T}
    finished::FinishedChannel{T}
    dag_worker::Task
    task_workers::Vector{Task}
    function PriorityScheduler{T}(sz = typemax(Int),background=false) where {T}
        if sz <= 0
            throw(ArgumentError("Scheduler buffer size must be a positive integer"))
        end
        dag            = DAG{T}(sz)
        nt             = Threads.nthreads()
        runnable       = RunnableChannel{T}()
        finished       = FinishedChannel{T}()
        dag_worker     = finished_to_runnable(dag,finished,runnable)
        task_workers   = consume_runnable(runnable,nt,background)
        sch            = new(dag,runnable,finished,dag_worker,task_workers)
        return sch
    end
end
PriorityScheduler(args...) = PriorityScheduler{DataFlowTask}(args...)

"""
    finished_to_runnable(dag,runnable,finished)

Worker which takes nodes from `finished`, remove them from the `dag`, and `put!`
new nodes in `runnable` if they become dependency-free.
"""
function finished_to_runnable(dag,finished,runnable)
    task = @async while true
        t    = take!(finished)
        # check if anyone in outlist is dependency-free after removal of t
        outlist = outneighbors(dag,t)
        remove_node!(dag,t)
        for j in outlist
            if isempty(inneighbors(dag,j))
                p = priority(j)
                put!(runnable,j,p)
            end
        end
    end
    errormonitor(task)
end

"""
    priority(t::DataFlowTask)

Function called to determine the scheduled priority of `t`. The default
imlementation simply retuns `t.priority`.
"""
function priority(j::DataFlowTask)
    j.priority
end

"""
    consume_runnable(runnable,nt)

Spawn `nt = Threads.nthreads()-1` background workers that will consume tasks
from `runnable` and execute them. The main thread (`Threads.threadid()==1`) is not used.
"""
function consume_runnable(runnable,nt,background)
    tasks = Task[]
    i1 = nt == 1 ? 1 : background ? 2 : 1
    for i in i1:nt
        t = @tspawnat i begin
            while true
                t    = take!(runnable)
                execute(t)
            end
        end
        push!(tasks,errormonitor(t))
    end
    return tasks
end

function execute(t::DataFlowTask)
    # FIXME: how to properly `execute` a Task without passing it to the
    # scheduler? Naive approach is to schedule and wait, but that seems silly in
    # the cases where where we know the task can just run immediately, so it is
    # safe to bypass the julia scheduler.
    schedule(t.task)
    wait(t)
    # task = t.task
    # task.result = Base.invokelatest(task.code)
    # task._state = 1
    # donenotify  = task.donenotify
    # lock(donenotify)
    # notify(donenotify)
    # unlock(donenotify)
end

# interface for PriorityScheduler
function spawn(tj,sch::PriorityScheduler)
    # each worker of the priority scheduler takes from runnable, so no need
    # to change sticky tj.task.sticky = false
    deps  = inneighbors(sch.dag,tj)
    if isempty(deps)
        p = priority(tj)
        put!(sch.runnable,tj,p)
    end
    # Since only tid 1 works on the dag, it is better to yield so that the
    # dag_worker can process finished nodes on the dag and make new nodes
    # runnable before continuing adding even more nodes to the dag
    yield()
    return tj
end

function Base.schedule(tj,sch::PriorityScheduler)
    deps  = inneighbors(sch.dag,tj)
    if isempty(deps)
        put!(sch.runnable,tj)
    end
    return tj
end


text:
using FaultTolerantControl
using Transducers
using Plots

function test()
    m = 1.0
    pos0 = zeros(3)
    g = 9.81
    controller = BacksteppingPositionControllerEnv(m)
    x0_controller = State(controller)(pos0, m, g)
    prob, df = sim(x0_controller,
                   apply_inputs(Dynamics!(controller),
                                pos_cmd=[2, 1, 3], Ṫd=1.0);
                   tf=10.0, savestep=0.01)
    ts = df.time
    xds = df.state |> Map(state -> state.ref_model.x_0) |> collect
    plot(ts, hcat(xds...)')
    # Tds = df.state |> Map(state -> state.Td) |> collect
    # plot(ts, hcat(Tds...)')
end


text:
# # Debug

Base.@deprecate debug(c) with_state_do(c)


text:
export AbstractInput,
    StartTargetInput,
    TargetSystemInput,
    HomotopyInput,
    ParameterSystemInput

const Inputs = Union{<:AbstractSystem, <:MPPolys, <:Composition}
const MPPolyInputs = Union{<:MPPolys, <:Composition}

const input_supported_keywords = [
    :parameters, :generic_parameters, :start_parameters, :target_parameters,
    :target_gamma, :start_gamma, :p₁, :p₀, :γ₁, :γ₀,
    # deprecated
    :startparameters, :targetparameters,
    :targetgamma, :startgamma,]

"""
    AbstractInput

An abstract type to model different input types.
"""
abstract type AbstractInput end

"""
    StartTargetInput(start::MPPolyInputs, target::MPPolyInputs, startsolutions)

Construct a `StartTargetInput` out of two systems of polynomials.
"""
struct StartTargetInput{P1<:Inputs, P2<:Inputs} <: AbstractInput
    start::P1
    target::P2

    function StartTargetInput{P1, P2}(start::P1, target::P2) where {P1<:Inputs, P2<:Inputs}
        if length(start) ≠ length(target)
            error("Cannot construct `StartTargetInput` since the lengths of `start` and `target` don't match.")
        end
        new(start, target)
    end
end
function StartTargetInput(start::P1, target::P2) where {P1<:Inputs, P2<:Inputs}
    StartTargetInput{P1, P2}(start, target)
end

"""
    HomotopyInput(H::Homotopy, startsolutions)

Construct a `HomotopyInput` for a homotopy `H` with given `startsolutions`.
"""
struct HomotopyInput{Hom<:AbstractHomotopy} <: AbstractInput
    H::Hom
end

"""
    TargetSystemInput(system::Inputs)

Construct a `TargetSystemInput`. This indicates that the system `system`
is the target system and a start system should be assembled.
"""
struct TargetSystemInput{S<:Inputs} <: AbstractInput
    system::S
end

"""
    ParameterSystemInput(F, parameters, p₁, p₀, startsolutions, γ₁=nothing, γ₂=nothing)

Construct a `ParameterSystemInput`.
"""
struct ParameterSystemInput{S<:Inputs} <: AbstractInput
    system::S
    parameters::Union{Nothing, Vector{<:MP.AbstractVariable}}
    p₁::AbstractVector
    p₀::AbstractVector
    γ₁::Union{Nothing, ComplexF64}
    γ₀::Union{Nothing, ComplexF64}
end

"""
    input_startsolutions(F::MPPolyInputs)
    input_startsolutions(F::AbstractSystem)
    input_startsolutions(G::MPPolyInputs, F::MPPolyInputs, startsolutions)
    input_startsolutions(F::MPPolyInputs, parameters, startsolutions; kwargs...)
    input_startsolutions(H::AbstractHomotopy, startsolutions)

Construct an `AbstractInput` and pass through startsolutions if provided.
Returns a named tuple `(input=..., startsolutions=...)`.
"""
function input_startsolutions(F::MPPolyInputs; parameters=nothing, kwargs...)
    # if parameters !== nothing this is actually the
    # input constructor for a parameter homotopy, but no startsolutions
    # are provided
    if parameters !== nothing
        startsolutions = [randn(ComplexF64, nvariables(F, parameters=parameters))]
        return input_startsolutions(F, startsolutions; parameters=parameters, kwargs...)
    end

    remove_zeros!(F)
    # square system and each polynomial is non-zero
    if length(F) == nvariables(F) && is_homogeneous(F)
        throw(ArgumentError("Cannot construct a start system for a square homogeneous system."))
    end
    if has_constant_polynomial(F)
        throw(ArgumentError("System contains a non-zero constant term"))
    end

    (input=TargetSystemInput(F), startsolutions=nothing)
end

function input_startsolutions(F::AbstractSystem)
    (input=TargetSystemInput(F), startsolutions=nothing)
end

function input_startsolutions(G::MPPolyInputs, F::MPPolyInputs, startsolutions=nothing)
    if length(G) ≠ length(F)
        throw(ArgumentError("Start and target system don't have the same length"))
    end
    check_zero_dimensional(F)
    if startsolutions === nothing
        startsolutions = [randn(ComplexF64, nvariables(F))]
    elseif isa(startsolutions, AbstractVector{<:Number})
        startsolutions = [startsolutions]
    end
    (input=StartTargetInput(G, F), startsolutions=startsolutions)
end

# need
input_startsolutions(F::MPPolyInputs, starts; kwargs...) = parameter_homotopy(F, starts; kwargs...)
input_startsolutions(F::AbstractSystem, starts; kwargs...) = parameter_homotopy(F, starts; kwargs...)

function parameter_homotopy(F::Inputs, startsolutions;
    parameters=(isa(F, AbstractSystem) ? nothing : error(ArgumentError("You need to pass `parameters=...` as a keyword argument."))),
    generic_parameters=nothing,
    start_parameters=generic_parameters, p₁ = start_parameters,
    target_parameters=generic_parameters, p₀ = target_parameters,
    start_gamma=nothing, γ₁ = start_gamma,
    target_gamma=nothing, γ₀ = target_gamma,
    # deprecated in 0.7
    startparameters=nothing,
    targetparameters=nothing,
    startgamma=nothing,
    targetgamma=nothing
    )

    # deprecation handling
    @deprecatekwarg startparameters start_parameters
    @deprecatekwarg targetparameters target_parameters
    @deprecatekwarg startgamma start_gamma
    @deprecatekwarg targetgamma target_gamma

    if γ₁ === nothing
        γ₁ = start_gamma
    end
    if γ₀ === nothing
        γ₀ = target_gamma
    end
    if p₁ === nothing
        p₁ = start_parameters
    end
    if p₀ === nothing
        p₀ = target_parameters
    end

    if p₁ === nothing
        error("You need to pass `generic_parameters=`, `start_parameters=` or `p₁=` as a keyword argument")
    elseif p₀ === nothing
        error("`target_parameters=` or `p₀=` need to be passed as a keyword argument.")
    end

    if length(p₁) != length(p₀) || (parameters !== nothing && length(parameters) != length(p₀))
        error("Number of parameters doesn't match!")
    end
    if startsolutions === nothing && parameters !== nothing
        startsolutions = [randn(ComplexF64, nvariables(F, parameters=parameters))]
    elseif isa(startsolutions, AbstractVector{<:Number})
        startsolutions = [startsolutions]
    end
    (input=ParameterSystemInput(F, parameters, p₁, p₀, γ₁, γ₀), startsolutions=startsolutions)
end

function input_startsolutions(H::AbstractHomotopy, startsolutions)
    if isa(startsolutions, AbstractVector{<:Number})
        startsolutions = [startsolutions]
    end
    (input=HomotopyInput(H), startsolutions=startsolutions)
end


text:
function solveSessionEvent(navAbilityClient::NavAbilityClient, client::Client)::String
    response = navAbilityClient.mutate(MutationOptions(
        "solveSession",
        MUTATION_SOLVESESSION,
        Dict(
            "client" => client,
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return "Error" end
    solveSession = get(data,"solveSession","Error")
    return solveSession
end

solveSession(navAbilityClient::NavAbilityClient, client::Client) = @async solveSessionEvent(navAbilityClient, client)

function solveFederatedEvent(navAbilityClient::NavAbilityClient, scope::Scope)::String
    response = navAbilityClient.mutate(MutationOptions(
        "solveFederated",
        MUTATION_SOLVEFEDERATED,
        Dict(
            "scope" => scope,
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return "Error" end
    solveSession = get(data,"solveSession","Error")
    return solveSession
end

solveFederated(navAbilityClient::NavAbilityClient, scope::Scope) = @async solveFederatedEvent(navAbilityClient, scope)

#


text:
# The following is largely inspired from JuMP/test/JuMPExtension.jl

@enum Level BOTH LOWER_ONLY UPPER_ONLY DUAL_OF_LOWER

abstract type AbstractBilevelModel <: JuMP.AbstractModel end

Base.broadcastable(model::AbstractBilevelModel) = Ref(model)

mutable struct BilevelModel <: AbstractBilevelModel
    # Structured data
    upper::JuMP.AbstractModel
    lower::JuMP.AbstractModel

    # Model data
    nextvaridx::Int                                 # Next variable index is nextvaridx+1
    variables::Dict{Int, JuMP.AbstractVariable}     # Map varidx -> variable
    varnames::Dict{Int, String}                     # Map varidx -> name
    varnames_rev::Dict{String, Int}                 # Map varidx -> name
    var_level::Dict{Int, Level}
    var_upper::Dict{Int, JuMP.AbstractVariableRef}
    var_lower::Dict{Int, JuMP.AbstractVariableRef} #
    var_info::Dict{Int, VariableInfo}

    # upper level decisions that are "parameters" of the second level
    upper_to_lower_link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef}
    # lower level decisions that are input to upper
    lower_to_upper_link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef}
    # lower level decisions that are input to upper
    upper_var_to_lower_ctr_link::Dict{JuMP.AbstractVariableRef, JuMP.ConstraintRef}
    # joint link
    link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef}

    nextconidx::Int                                 # Next constraint index is nextconidx+1
    constraints::Dict{Int, JuMP.AbstractConstraint} # Map conidx -> variable
    connames::Dict{Int, String}                     # Map varidx -> name
    connames_rev::Dict{String, Int}                     # Map varidx -> name
    ctr_level::Dict{Int, Level}
    ctr_upper::Dict{Int, JuMP.ConstraintRef}
    ctr_lower::Dict{Int, JuMP.ConstraintRef}
    ctr_info::Dict{Int, ConstraintInfo}

    upper_objective_sense::MOI.OptimizationSense
    upper_objective_function::JuMP.AbstractJuMPScalar

    lower_objective_sense::MOI.OptimizationSense
    lower_objective_function::JuMP.AbstractJuMPScalar

    # solution data
    need_rebuild::Bool
    need_rebuild_names_ctr::Bool
    need_rebuild_names_var::Bool
    solver#::MOI.ModelLike
    mode
    upper_to_sblm
    lower_to_sblm
    lower_dual_to_sblm
    sblm_to_solver
    lower_primal_dual_map

    solve_time::Float64
    build_time::Float64

    copy_names::Bool
    copy_names_to_solver::Bool
    pass_start::Bool

    objdict::Dict{Symbol, Any}    # Same that JuMP.Model's field `objdict`

    function BilevelModel()

        model = new(
            JuMP.Model(),
            JuMP.Model(),

            # var
            0, Dict{Int, JuMP.AbstractVariable}(),
            Dict{Int, String}(),Dict{String, Int}(),
            Dict{Int, Level}(), Dict{Int, JuMP.AbstractVariable}(), Dict{Int, JuMP.AbstractVariable}(),
            Dict{Int, VariableInfo}(),
            # links
            Dict{JuMP.AbstractVariable, JuMP.AbstractVariable}(), Dict{JuMP.AbstractVariable, JuMP.AbstractVariable}(),
            Dict{JuMP.AbstractVariable, JuMP.ConstraintRef}(),
            Dict{JuMP.AbstractVariable, JuMP.AbstractVariable}(),
            #ctr
            0, Dict{Int, JuMP.AbstractConstraint}(),
            Dict{Int, String}(),Dict{String, Int}(),
            Dict{Int, Level}(), Dict{Int, JuMP.AbstractConstraint}(), Dict{Int, JuMP.AbstractConstraint}(),
            Dict{Int, ConstraintInfo}(),
            #obj
            MOI.FEASIBILITY_SENSE,
            zero(JuMP.GenericAffExpr{Float64, BilevelVariableRef}), # Model objective
            MOI.FEASIBILITY_SENSE,
            zero(JuMP.GenericAffExpr{Float64, BilevelVariableRef}), # Model objective

            true,
            true,
            true,
            nothing,
            NoMode{Float64},
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,

            NaN,
            NaN,
            false,
            false,
            true,
            Dict{Symbol, Any}(),
            )

        return model
    end
end
function BilevelModel(optimizer_constructor;
    mode::AbstractBilevelSolverMode = SOS1Mode(),
    bridge_constraints::Bool=true,
    with_names=true)
    bm = BilevelModel()
    bm.copy_names=with_names
    set_mode(bm, mode)
    JuMP.set_optimizer(bm, optimizer_constructor;
        bridge_constraints=bridge_constraints,
        with_names=with_names)
    return bm
end
function set_mode(bm::BilevelModel, mode::AbstractBilevelSolverMode)
    bm.mode = deepcopy(mode)
    reset!(bm.mode)
    return bm
end

abstract type InnerBilevelModel <: AbstractBilevelModel end
struct UpperModel <: InnerBilevelModel
    m::BilevelModel
end
Upper(m::BilevelModel) = UpperModel(m)
UpperToLower(m::BilevelModel) = UpperModel(m)
struct LowerModel <: InnerBilevelModel
    m::BilevelModel
end
Lower(m::BilevelModel) = LowerModel(m)
LowerToUpper(m::BilevelModel) = LowerModel(m)
bilevel_model(m::InnerBilevelModel) = m.m
mylevel_model(m::UpperModel) = bilevel_model(m).upper
mylevel_model(m::LowerModel) = bilevel_model(m).lower
level(m::LowerModel) = LOWER_ONLY
level(m::UpperModel) = UPPER_ONLY
mylevel_ctr_list(m::LowerModel) = bilevel_model(m).ctr_lower
mylevel_ctr_list(m::UpperModel) = bilevel_model(m).ctr_upper
mylevel_var_list(m::LowerModel) = bilevel_model(m).var_lower
mylevel_var_list(m::UpperModel) = bilevel_model(m).var_upper

# obj

mylevel_obj_sense(m::LowerModel) = bilevel_model(m).lower_objective_sense
mylevel_obj_function(m::LowerModel) = bilevel_model(m).lower_objective_function
mylevel_obj_sense(m::UpperModel) = bilevel_model(m).upper_objective_sense
mylevel_obj_function(m::UpperModel) = bilevel_model(m).upper_objective_function

set_mylevel_obj_sense(m::LowerModel, val) = bilevel_model(m).lower_objective_sense = val
set_mylevel_obj_function(m::LowerModel, val) = bilevel_model(m).lower_objective_function = val
set_mylevel_obj_sense(m::UpperModel, val) = bilevel_model(m).upper_objective_sense = val
set_mylevel_obj_function(m::UpperModel, val) = bilevel_model(m).upper_objective_function = val

function set_link!(m::UpperModel, upper::JuMP.AbstractVariableRef, lower::JuMP.AbstractVariableRef)
    bilevel_model(m).upper_to_lower_link[upper] = lower
    bilevel_model(m).link[upper] = lower
    nothing
end
function set_link!(m::LowerModel, upper::JuMP.AbstractVariableRef, lower::JuMP.AbstractVariableRef)
    bilevel_model(m).lower_to_upper_link[lower] = upper
    bilevel_model(m).link[upper] = lower
    nothing
end

# Models to deal with variables that are not exchanged between models
abstract type SingleBilevelModel <: AbstractBilevelModel end
struct UpperOnlyModel <: SingleBilevelModel
    m::BilevelModel
end
UpperOnly(m::BilevelModel) = UpperOnlyModel(m)
struct LowerOnlyModel <: SingleBilevelModel
    m::BilevelModel
end
LowerOnly(m::BilevelModel) = LowerOnlyModel(m)

bilevel_model(m::SingleBilevelModel) = m.m
mylevel_model(m::UpperOnlyModel) = bilevel_model(m).upper
mylevel_model(m::LowerOnlyModel) = bilevel_model(m).lower
level(::LowerOnlyModel) = LOWER_ONLY
level(::UpperOnlyModel) = UPPER_ONLY
mylevel_var_list(m::LowerOnlyModel) = bilevel_model(m).var_lower
mylevel_var_list(m::UpperOnlyModel) = bilevel_model(m).var_upper

in_upper(l::Level) = l == BOTH || l == UPPER_ONLY || l == DUAL_OF_LOWER
in_lower(l::Level) = l == BOTH || l == LOWER_ONLY

#### Model ####

# Variables
struct BilevelVariableRef <: JuMP.AbstractVariableRef
    model::BilevelModel # `model` owning the variable
    idx::Int       # Index in `model.variables`
    level::Level
end
function BilevelVariableRef(model, idx)
    return BilevelVariableRef(model, idx, model.var_level[idx])
end

# Constraints
const BilevelConstraintRef = JuMP.ConstraintRef{BilevelModel, Int}#, Shape <: AbstractShape

# Objective

# Etc

JuMP.object_dictionary(m::BilevelModel) = m.objdict
JuMP.object_dictionary(m::AbstractBilevelModel) = JuMP.object_dictionary(bilevel_model(m))

function JuMP.index(d::Dict)
    ret = Dict{VI,VI}()
    # sizehint!(ret, length(d))
    for (k,v) in d
        ret[JuMP.index(k)] = JuMP.index(v)
    end
    return ret
end
function index2(d::Dict)
    ret = Dict{VI,CI}()
    # sizehint!(ret, length(d))
    for (k,v) in d
        ret[JuMP.index(k)] = JuMP.index(v)
    end
    return ret
end

# Names
JuMP.name(vref::BilevelVariableRef) = vref.model.varnames[vref.idx]
function JuMP.set_name(vref::BilevelVariableRef, name::String)
    vref.model.varnames[vref.idx] = name
end
JuMP.name(cref::BilevelConstraintRef) = cref.model.connames[cref.index]
function JuMP.set_name(cref::BilevelConstraintRef, name::String)
    cref.model.connames[cref.index] = name
end
function JuMP.variable_by_name(model::BilevelModel, name::String)
    if model.need_rebuild_names_var
        model.varnames_rev = Dict(v => k for (k,v) in model.varnames)
    end
    idx = model.varnames_rev[name]
    return BilevelVariableRef(model, idx)
end
function JuMP.constraint_by_name(model::BilevelModel, name::String)
    if model.need_rebuild_names_ctr
        model.connames_rev = Dict(v => k for (k,v) in model.connames)
    end
    idx = model.connames_rev[name]
    return BilevelConstraintRef(model, idx)
end

function JuMP.primal_status(model::BilevelModel)
    _check_solver(model)
    return MOI.get(model.solver, MOI.PrimalStatus())
end
function JuMP.primal_status(model::InnerBilevelModel)
    return JuMP.primal_status(model.m)
end

# Statuses

JuMP.dual_status(::BilevelModel) = error(
    "Dual status cant be queried for BilevelModel, but you can query for Upper and Lower models.")
function JuMP.dual_status(model::UpperModel)
    _check_solver(model.m)
    return MOI.get(model.m.solver, MOI.DualStatus())
end
function JuMP.dual_status(model::LowerModel)
    _check_solver(model.m)
    return MOI.get(model.m.solver, MOI.PrimalStatus())
end

function JuMP.termination_status(model::BilevelModel)
    _check_solver(model)
    return MOI.get(model.solver, MOI.TerminationStatus())
end
function JuMP.raw_status(model::BilevelModel)
    _check_solver(model)
    return MOI.get(model.solver, MOI.RawStatusString())
end

# Replace variables

function replace_variables(var::BilevelVariableRef,
    model::BilevelModel,
    inner::JuMP.AbstractModel,
    variable_map::Dict{Int, V},
    level::Level) where {V<:JuMP.AbstractVariableRef}
    if var.model === model && in_level(var, level)
        return variable_map[var.idx]
    elseif var.model === model
        error("Variable $(var) belonging Only to $(var.level) level, was added in the $(level) level.")
    else
        error("A BilevelModel cannot have expression using variables of a BilevelModel different from itself")
    end
end
function replace_variables(aff::JuMP.GenericAffExpr{C, BilevelVariableRef},
    model::BilevelModel,
    inner::JuMP.AbstractModel,
    variable_map::Dict{Int, V},
    level::Level) where {C,V<:JuMP.AbstractVariableRef}
    result = JuMP.GenericAffExpr{C, JuMP.VariableRef}(0.0)#zero(aff)
    result.constant = aff.constant
    for (coef, var) in JuMP.linear_terms(aff)
        JuMP.add_to_expression!(result,
        coef,
        replace_variables(var, model, model, variable_map, level))
    end
    return result
end
function replace_variables(quad::JuMP.GenericQuadExpr{C, BilevelVariableRef},
    model::BilevelModel,
    inner::JuMP.AbstractModel,
    variable_map::Dict{Int, V},
    level::Level) where {C,V<:JuMP.AbstractVariableRef}
    aff = replace_variables(quad.aff, model, model, variable_map, level)
    quadv = JuMP.GenericQuadExpr{C, JuMP.VariableRef}(aff)
    for (coef, var1, var2) in JuMP.quad_terms(quad)
        JuMP.add_to_expression!(quadv,
        coef,
        replace_variables(var1, model, model, variable_map, level),
        replace_variables(var2, model, model, variable_map, level))
    end
    return quadv
end
replace_variables(funcs::Vector, args...) = map(f -> replace_variables(f, args...), funcs)

function print_lp(m, name, file_format = MOI.FileFormats.FORMAT_AUTOMATIC)
    dest = MOI.FileFormats.Model(format = file_format, filename = name)
    MOI.copy_to(dest, m)
    MOI.write_to_file(dest, name)
end

# Optimize

JuMP.optimize!(::T) where {T<:AbstractBilevelModel} =
    error("Can't solve a model of type: $T ")
function JuMP.optimize!(model::BilevelModel;
    lower_prob = "", upper_prob = "", bilevel_prob = "", solver_prob = "", file_format = MOI.FileFormats.FORMAT_AUTOMATIC)

    if model.mode === nothing
        error("No solution mode selected, use `set_mode(model, mode)` or initialize with `BilevelModel(optimizer_constructor, mode = some_mode)`")
    else
        mode = model.mode
    end

    _check_solver(model)

    solver = model.solver #optimizer#MOI.Bridges.full_bridge_optimizer(optimizer, Float64)
    if true#!MOI.is_empty(solver)
        MOI.empty!(solver)
    end

    upper = JuMP.backend(model.upper)
    lower = JuMP.backend(model.lower)

    if length(lower_prob) > 0
        print_lp(lower, lower_prob, file_format)
    end
    if length(upper_prob) > 0
        print_lp(upper, upper_prob, file_format)
    end

    t0 = time()

    moi_upper = JuMP.index.(
        collect(values(model.upper_to_lower_link)))
    moi_link = JuMP.index(model.link)
    moi_link2 = index2(model.upper_var_to_lower_ctr_link)

    # build bound for FortunyAmatMcCarlMode
    build_bounds!(model, mode)

    single_blm, upper_to_sblm, lower_to_sblm, lower_primal_dual_map, lower_dual_to_sblm =
        build_bilevel(upper, lower, moi_link, moi_upper, mode, moi_link2,
            copy_names = model.copy_names, pass_start = model.pass_start)

    # pass lower level dual variables info (start, upper, lower)
    for (idx, info) in model.ctr_info
        if haskey(model.ctr_lower, idx)
            ctr = model.ctr_lower[idx]
            pre_duals = lower_primal_dual_map.primal_con_dual_var[JuMP.index(ctr)] # vector
            duals = map(x->lower_dual_to_sblm[x], pre_duals)
            pass_dual_info(single_blm, duals, info)
        end
    end
    # pass lower & upper level primal variables info (upper, lower)
    for (idx, info) in model.var_info
        if haskey(model.var_lower, idx)
            var = lower_to_sblm[JuMP.index(model.var_lower[idx])]
        elseif haskey(model.var_upper, idx)
            var = upper_to_sblm[JuMP.index(model.var_upper[idx])]
        else
            continue
        end
        pass_primal_info(single_blm, var, info)
    end
    if length(bilevel_prob) > 0
        print_lp(single_blm, bilevel_prob, file_format)
    end

    sblm_to_solver = MOI.copy_to(solver, single_blm, copy_names = model.copy_names_to_solver)

    if length(solver_prob) > 0
        print_lp(solver, solver_prob, file_format)
    end

    model.upper_to_sblm = upper_to_sblm
    model.lower_to_sblm = lower_to_sblm
    model.lower_dual_to_sblm = lower_dual_to_sblm
    model.lower_primal_dual_map = lower_primal_dual_map
    model.sblm_to_solver = sblm_to_solver

    t1 = time()
    model.build_time = t1 - t0

    MOI.optimize!(solver)

    model.solve_time = time() - t1

    return nothing
end

# Extra info

function pass_primal_info(single_blm, primal, info::VariableInfo)
    if !isnan(info.upper) &&
        !MOI.is_valid(single_blm, CI{SVF,LT{Float64}}(primal.value))
        MOI.add_constraint(single_blm,
            SVF(primal), LT{Float64}(info.upper))
    end
    if !isnan(info.lower) &&
        !MOI.is_valid(single_blm, CI{SVF,GT{Float64}}(primal.value))
        MOI.add_constraint(single_blm,
            SVF(primal), GT{Float64}(info.lower))
    end
    return
end

function pass_dual_info(single_blm, dual, info::ConstraintInfo{Float64})
    if !isnan(info.start)
        MOI.set(single_blm, MOI.VariablePrimalStart(), dual[], info.start)
    end
    if !isnan(info.upper) &&
        !MOI.is_valid(single_blm, CI{SVF,LT{Float64}}(dual[].value))
        MOI.add_constraint(single_blm,
            SVF(dual[]), LT{Float64}(info.upper))
    end
    if !isnan(info.lower) &&
        !MOI.is_valid(single_blm, CI{SVF,GT{Float64}}(dual[].value))
        MOI.add_constraint(single_blm,
            SVF(dual[]), GT{Float64}(info.lower))
    end
    return
end
function pass_dual_info(single_blm, dual, info::ConstraintInfo{Vector{Float64}})
    for i in eachindex(dual)
        if !isnan(info.start[i])
            MOI.set(single_blm, MOI.VariablePrimalStart(), dual[i], info.start[i])
        end
        if !isnan(info.upper[i]) &&
            !MOI.is_valid(single_blm, CI{SVF,LT{Float64}}(dual[i].value))
            MOI.add_constraint(single_blm,
                SVF(dual[i]), LT{Float64}(info.upper[i]))
        end
        if !isnan(info.lower[i]) &&
            !MOI.is_valid(single_blm, CI{SVF,GT{Float64}}(dual[i].value))
            MOI.add_constraint(single_blm,
                SVF(dual[i]), MOI.GreaterThan{Float64}(info.lower[i]))
        end
    end
    return
end

# Bounds

function build_bounds!(::BilevelModel, ::AbstractBilevelSolverMode{T}) where T
    return nothing
end
function build_bounds!(model::BilevelModel, mode::AbstractBoundedMode{T}) where T
    return _build_bounds!(model, mode.cache)
end
function _build_bounds!(model::BilevelModel, mode::ComplementBoundCache)
    # compute variable bounds for FA mode
    fa_vi_up = mode.upper
    fa_vi_lo = mode.lower
    fa_vi_ld = mode.ldual
    empty!(fa_vi_up)
    empty!(fa_vi_lo)
    empty!(fa_vi_ld)
    for (idx, _info) in model.var_info
        if haskey(model.var_lower, idx)
            var = JuMP.index(model.var_lower[idx])
            is_lower = true
        elseif haskey(model.var_upper, idx)
            var = JuMP.index(model.var_upper[idx])
            is_lower = false
        else
            continue
        end
        vref = BilevelVariableRef(model, idx)

        is_dual = vref.level == DUAL_OF_LOWER

        info = deepcopy(_info)

        lb = JuMP.has_lower_bound(vref) ? JuMP.lower_bound(vref) : -Inf
        ub = JuMP.has_upper_bound(vref) ? JuMP.upper_bound(vref) : +Inf
        info.upper = min(ub, inf_if_nan(+, info.upper))
        info.lower = max(lb, inf_if_nan(-, info.lower))
        if is_lower
            fa_vi_lo[var] = info
        else
            fa_vi_up[var] = info
        end
        @assert info.lower <= info.upper
    end
    for (idx, _info) in model.ctr_info
        if haskey(model.ctr_lower, idx)
            ctr = JuMP.index(model.ctr_lower[idx])
            info = deepcopy(_info)
            cref = BilevelConstraintRef(model, idx)
            # scalar
            ub = dual_upper_bound(ctr)
            lb = dual_lower_bound(ctr)
            info.upper = min.(ub, inf_if_nan.(+, info.upper))
            info.lower = max.(lb, inf_if_nan.(-, info.lower))
            @assert sum(info.lower .<= info.upper) == length(info.upper)
            # TODO vector
            fa_vi_ld[ctr] = info
        end
    end
    return nothing
end

inf_if_nan(::typeof(+), val) = ifelse(isnan(val), Inf, val)
inf_if_nan(::typeof(-), val) = ifelse(isnan(val), -Inf, val)

dual_lower_bound(::CI{F,LT{T}}) where {F,T} = -Inf
dual_upper_bound(::CI{F,LT{T}}) where {F,T} =  0.0

dual_lower_bound(::CI{F,GT{T}}) where {F,T} =  0.0
dual_upper_bound(::CI{F,GT{T}}) where {F,T} = +Inf

dual_lower_bound(::CI{F,ET{T}}) where {F,T} =  0.0
dual_upper_bound(::CI{F,ET{T}}) where {F,T} =  0.0

dual_lower_bound(::CI{F,S}) where {F,S} = -Inf
dual_upper_bound(::CI{F,S}) where {F,S} = +Inf

# Initialize

function _check_solver(bm::BilevelModel)
    if bm.solver === nothing
        error("No solver attached, use `set_optimizer(model, optimizer_constructor)` or initialize with `BilevelModel(optimizer_constructor)`")
    end
end

function JuMP.set_optimizer(bm::BilevelModel, optimizer_constructor;
    bridge_constraints::Bool=true, with_names=true)
    # error_if_direct_mode(model, :set_optimizer)
    if bridge_constraints
        # We set `with_names=false` because the names are handled by the first
        # caching optimizer. If `default_copy_to` without names is supported,
        # no need for a second cache.
        optimizer = MOI.instantiate(optimizer_constructor, with_bridge_type=Float64, with_names=with_names)
        # for bridge_type in model.bridge_types
        #     _moi_add_bridge(optimizer, bridge_type)
        # end
    else
        optimizer = MOI.instantiate(optimizer_constructor)
    end

    bm.solver = optimizer
    if !MOI.is_empty(bm.solver)
        error("Calling the `optimizer_constructor` must return an empty optimizer")
    end
    return bm
end

function pass_cache(bm::BilevelModel, mode::FortunyAmatMcCarlMode{T}) where T
    mode.cache = bm.mode.cache
    return nothing
end
function pass_cache(bm::BilevelModel, mode::AbstractBilevelSolverMode{T}) where T
    return nothing
end

function check_mixed_mode(::MixedMode{T}) where T
end
function check_mixed_mode(mode)
    error("Cant set/get mode on a specific object because the base mode is $mode while it should be MixedMode in this case. Run `set_mode(model, BilevelJuMP.MixedMode())`")
end
function set_mode(ci::BilevelConstraintRef, mode::AbstractBilevelSolverMode{T}) where T
    bm = ci.model
    check_mixed_mode(bm.mode)
    _mode = deepcopy(mode)
    pass_cache(bm, _mode)
    ctr = JuMP.index(bm.ctr_lower[ci.index])
    bm.mode.constraint_mode_map_c[ctr] = _mode
    return nothing
end
function unset_mode(ci::BilevelConstraintRef)
    bm = ci.model
    check_mixed_mode(bm.mode)
    ctr = JuMP.index(bm.ctr_lower[ci.index])
    delete!(bm.mode.constraint_mode_map_c, ctr)
    return nothing
end

function get_mode(ci::BilevelConstraintRef)
    bm = ci.model
    check_mixed_mode(bm.mode)
    ctr = JuMP.index(bm.ctr_lower[ci.index])
    if haskey(bm.mode.constraint_mode_map_c, ctr)
        return bm.mode.constraint_mode_map_c[ctr]
    else
        return nothing
    end
end

function set_mode(::BilevelConstraintRef, ::MixedMode{T}) where T
    error("Cant set MixedMode in a specific constraint")
end
function set_mode(::BilevelConstraintRef, ::StrongDualityMode{T}) where T
    error("Cant set StrongDualityMode in a specific constraint")
end

function set_mode(vi::BilevelVariableRef, mode::AbstractBilevelSolverMode{T}) where T
    bm = vi.model
    check_mixed_mode(bm.mode)
    _mode = deepcopy(mode)
    pass_cache(bm, _mode)
    var = JuMP.index(bm.var_lower[vi.idx])
    # var is the MOI backend index in jump
    bm.mode.constraint_mode_map_v[var] = _mode
    return nothing
end
function unset_mode(vi::BilevelVariableRef)
    bm = vi.model
    check_mixed_mode(bm.mode)
    var = JuMP.index(bm.var_lower[vi.idx])
    delete!(bm.mode.constraint_mode_map_v, var)
    return nothing
end

function get_mode(vi::BilevelVariableRef)
    bm = vi.model
    check_mixed_mode(bm.mode)
    var = JuMP.index(bm.var_lower[vi.idx])
    if haskey(bm.mode.constraint_mode_map_v, var)
        return bm.mode.constraint_mode_map_v[var]
    else
        return nothing
    end
end

function set_mode(::BilevelVariableRef, ::MixedMode{T}) where T
    error("Cant set MixedMode in a specific variable")
end
function set_mode(::BilevelVariableRef, ::StrongDualityMode{T}) where T
    error("Cant set StrongDualityMode in a specific variable")
end


text:
export PointNeuron, WhiteSteady, WhiteCmod, WhiteVmod

@with_kw struct PointNeuron
  τ_v::Float64 # membrane time constant, ms
end

abstract type DriveParams end

@with_kw struct WhiteSteady <: DriveParams
  μ::Float64
  σ_w::Float64
end

function white_drive(I_args::WhiteSteady, dt, t)
  @unpack μ, σ_w = I_args
  μ*dt+σ_w*sqrt(2*dt)*randn()
end

@with_kw struct WhiteCmod <: DriveParams
  μ::Float64
  σ_w::Float64
  ϵ_c::Float64
  Ω::Float64
end

function white_drive(I_args::WhiteCmod, dt, t)
  @unpack μ,σ_w, ϵ_c, Ω = I_args
  μ*dt+σ_w*sqrt(2*dt)*randn()+ϵ_c*dt*sin(Ω*t)
end

@with_kw struct WhiteVmod <: DriveParams
  μ::Float64
  σ_w::Float64
  ϵ_v::Float64
  Ω::Float64
end

function white_drive(I_args::WhiteVmod, dt, t)
  @unpack μ, σ_w, ϵ_v, Ω = I_args
  μ*dt+σ_w*(1+ϵ_v*sin(Ω*t))*sqrt(2*dt)*randn()
end

function point_white_vtime(N_args::PointNeuron, I_args::DriveParams, T::TimeAxis; seed = 0)
	if seed != 0
	  Random.seed!(seed)
	end
  @unpack τ_v = N_args
  @unpack dt, N = T

	V = zeros(N)
  for i=2:N
  	V[i] = V[i-1]*(1-dt)+white_drive(I_args,dt,(i-1)*dt)
  end
  t = range(0,step=dt,length=N)
  return t.*τ_v, V
end
export point_white_vtime

function point_white_steady_test()
  T = TimeAxis(dt = 0.05, N = 1000)
  example = PointNeuron(τ_v = 10.0)
  drive = WhiteSteady(μ=5.0, σ_w = 1.0)
  @time t, V = point_white_vtime(example, drive, T; seed = 1000)
  return t, V
end
export point_white_steady_test

function point_white_cmod_test()
  T = TimeAxis(dt = 0.05, N = 1000)
  example = PointNeuron(τ_v = 10.0)
  drive = WhiteCmod(μ=5.0, σ_w = 1.0, ϵ_c = 5.0, Ω = 0.1*pi)
  @time t, V = point_white_vtime(example, drive, T; seed = 1000)
  return t, V
end
export point_white_cmod_test

function point_white_vmod_test()
  T = TimeAxis(dt = 0.05, N = 1000)
  example = PointNeuron(τ_v = 10.0)
  drive = WhiteVmod(μ=5.0, σ_w = 1.0, ϵ_v = 5.0, Ω = 0.1*pi)
  @time t, V = point_white_vtime(example, drive, T; seed = 1000)
  return t, V
end
export point_white_vmod_test

max_stars_repo_name: rpgowers/NoisyNeuron

max_issues_repo_name: rpgowers/NoisyNeuron

max_forks_repo_name: rpgowers/NoisyNeuron