text:
# MIP01 from CPLEX.jl
function int1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # an example on mixed integer programming
    #
    #   maximize 1.1x + 2 y + 5 z
    #
    #   s.t.  x + y + z <= 10
    #         x + 2 y + z <= 15
    #
    #         x is continuous: 0 <= x <= 5
    #         y is integer: 0 <= y <= 10
    #         z is binary

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Integer)

    MOI.empty!(model)
    @test MOI.is_empty(model)

    v = MOI.add_variables(model, 3)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3

    cf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0,1.0], v), 0.0)
    c = MOI.add_constraint(model, cf, MOI.LessThan(10.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,2.0,1.0], v), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.LessThan(15.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 2

    vc1 = MOI.add_constraint(model, MOI.SingleVariable(v[1]), MOI.Interval(0.0, 5.0))
    # We test this after the creation of every `SingleVariable` constraint
    # to ensure a good coverage of corner cases.
    @test vc1.value == v[1].value
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Interval{Float64}}()) == 1

    vc2 = MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.Interval(0.0, 10.0))
    @test vc2.value == v[2].value
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Interval{Float64}}()) == 2
    vc3 = MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.Integer())
    @test vc3.value == v[2].value
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Integer}()) == 1

    vc4 = MOI.add_constraint(model, MOI.SingleVariable(v[3]), MOI.ZeroOne())
    @test vc4.value == v[3].value
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.ZeroOne}()) == 1

    objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.1, 2.0, 5.0], v), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) in [ MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT ]

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 19.4 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4,5,1] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 10 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c2) ≈ 15 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) >= 19.4 - atol
        # FIXME the following are currently not implemented in MockOptimizer
#        @test MOI.get(model, MOI.RelativeGap()) >= 0.0
#        @test MOI.get(model, MOI.SolveTime()) >= 0.0
#        @test MOI.get(model, MOI.SimplexIterations()) >= 0
#        @test MOI.get(model, MOI.BarrierIterations()) >= 0
#        @test MOI.get(model, MOI.NodeCount()) >= 0
    end
end

# sos from CPLEX.jl" begin
function int2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "SOSI" begin
        @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
        @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test MOI.supports(model, MOI.ObjectiveSense())
        @test MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.SOS1{Float64})
        @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.LessThan{Float64})

        MOI.empty!(model)
        @test MOI.is_empty(model)

        v = MOI.add_variables(model, 3)
        @test MOI.get(model, MOI.NumberOfVariables()) == 3
        vc1 = MOI.add_constraint(model, MOI.SingleVariable(v[1]), MOI.LessThan(1.0))
        @test vc1.value == v[1].value
        vc2 = MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.LessThan(1.0))
        @test vc2.value == v[2].value
        vc3 = MOI.add_constraint(model, MOI.SingleVariable(v[3]), MOI.LessThan(2.0))
        @test vc3.value == v[3].value

        c1 = MOI.add_constraint(model, MOI.VectorOfVariables([v[1], v[2]]), MOI.SOS1([1.0, 2.0]))
        c2 = MOI.add_constraint(model, MOI.VectorOfVariables([v[1], v[3]]), MOI.SOS1([1.0, 2.0]))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SOS1{Float64}}()) == 2

        #=
            To allow for permutations in the sets and variable vectors
            we're going to sort according to the weights
        =#
        cs_sos = MOI.get(model, MOI.ConstraintSet(), c2)
        cf_sos = MOI.get(model, MOI.ConstraintFunction(), c2)
        p = sortperm(cs_sos.weights)
        @test cs_sos.weights[p] ≈ [1.0, 2.0] atol=atol rtol=rtol
        @test cf_sos.variables[p] == v[[1,3]]

        objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0, 1.0], v), 0.0)
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

        if config.solve
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

            MOI.optimize!(model)

            @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3 atol=atol rtol=rtol

            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0,1,2] atol=atol rtol=rtol
        end

        MOI.delete(model, c1)
        MOI.delete(model, c2)

        if config.solve
            MOI.optimize!(model)

            @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 5 atol=atol rtol=rtol

            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1,1,2] atol=atol rtol=rtol
        end
    end
    @testset "SOSII" begin
        @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
        @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test MOI.supports(model, MOI.ObjectiveSense())
        @test MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.SOS1{Float64})
        @test MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.SOS2{Float64})
        @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.ZeroOne)
        @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})

        MOI.empty!(model)
        @test MOI.is_empty(model)

        v = MOI.add_variables(model, 10)
        @test MOI.get(model, MOI.NumberOfVariables()) == 10

        bin_constraints = []
        for i in 1:8
            vc = MOI.add_constraint(model, MOI.SingleVariable(v[i]), MOI.Interval(0.0, 2.0))
            @test vc.value == v[i].value
            push!(bin_constraints, MOI.add_constraint(model, MOI.SingleVariable(v[i]), MOI.ZeroOne()))
            @test bin_constraints[i].value == v[i].value
        end

        MOI.add_constraint(model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,2.0,3.0,-1.0], v[[1,2,3,9]]), 0.0),
            MOI.EqualTo(0.0)
        )

        MOI.add_constraint(model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5.0,4.0,7.0,2.0,1.0,-1.0], v[[4,5,6,7,8,10]]), 0.0),
            MOI.EqualTo(0.0)
        )

        MOI.add_constraint(model,
            MOI.VectorOfVariables(v[[1, 2, 3]]),
            MOI.SOS1([1.0, 2.0, 3.0])
        )

        vv   = MOI.VectorOfVariables(v[[4,5,6,7,8]])
        sos2 = MOI.SOS2([5.0, 4.0, 7.0, 2.0, 1.0])
        c = MOI.add_constraint(model, vv, sos2)

        #=
            To allow for permutations in the sets and variable vectors
            we're going to sort according to the weights
        =#
        cs_sos = MOI.get(model, MOI.ConstraintSet(), c)
        cf_sos = MOI.get(model, MOI.ConstraintFunction(), c)
        p = sortperm(cs_sos.weights)
        @test cs_sos.weights[p] ≈ [1.0, 2.0, 4.0, 5.0, 7.0] atol=atol rtol=rtol
        @test cf_sos.variables[p] == v[[8,7,5,4,6]]

        objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [v[9], v[10]]), 0.0)
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

        if config.solve
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

            MOI.optimize!(model)

            @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 15.0 atol=atol rtol=rtol

            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 3.0, 12.0] atol=atol rtol=rtol
        end

        for cref in bin_constraints
            MOI.delete(model, cref)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 30.0 atol=atol rtol=rtol

            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0.0, 0.0, 2.0, 2.0, 0.0, 2.0, 0.0, 0.0, 6.0, 24.0] atol=atol rtol=rtol
        end
    end
end

# CPLEX #76
function int3test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # integer knapsack problem
    # max   z - 0.5 ( b1 + b2 + b3) / 40
    # s.t.  0 <= z - 0.5 eᵀ b / 40 <= 0.999
    #       b1, b2, ... b10 ∈ {0, 1}
    #       z in {0, 1, 2, ..., 100}

    MOI.empty!(model)
    @test MOI.is_empty(model)

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})

    z = MOI.add_variable(model)
    vc1 = MOI.add_constraint(model, MOI.SingleVariable(z), MOI.Integer())
    @test vc1.value == z.value
    vc2 = MOI.add_constraint(model, MOI.SingleVariable(z), MOI.Interval(0.0, 100.0))
    @test vc2.value == z.value

    b = MOI.add_variables(model, 10)

    for bi in b
        vc = MOI.add_constraint(model, MOI.SingleVariable(bi), MOI.ZeroOne())
        @test vc.value == bi.value
    end

    c = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(vcat(1.0, fill(-0.5 / 40, 10)), vcat(z, b)), 0.0), MOI.Interval(0.0, 0.999))

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(vcat(1.0, fill(-0.5 / 40, 3)), vcat(z, b[1:3])), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1 atol=atol rtol=rtol

        # test for CPLEX.jl #76
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1 atol=atol rtol=rtol
    end
end

# Mixed-integer linear problems

function knapsacktest(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # integer knapsack problem
    # max 5a + 3b + 2c + 7d + 4e
    # st  2a + 8b + 4c + 2d + 5e <= 10
    #                  a,b,c,d,e ∈ binary

    MOI.empty!(model)
    @test MOI.is_empty(model)

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})

    v = MOI.add_variables(model, 5)
    @test MOI.get(model, MOI.NumberOfVariables()) == 5

    for vi in v
        vc = MOI.add_constraint(model, MOI.SingleVariable(vi), MOI.ZeroOne())
        @test vc.value == vi.value
    end
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.ZeroOne}()) == 5
    c = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 8.0, 4.0, 2.0, 5.0], v), 0.0), MOI.LessThan(10.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5.0, 3.0, 2.0, 7.0, 4.0], v), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    if MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
        MOI.set(model, MOI.VariablePrimalStart(), v, [0.0, 0.0, 0.0, 0.0, 0.0])
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) in [ MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT ]

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 16 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1, 0, 0, 1, 1] atol=atol rtol=rtol
    end
end

function indicator1_test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # linear problem with indicator constraint
    # max  2x1 + 3x2
    # s.t. x1 + x2 <= 10
    #      z1 ==> x2 <= 8
    #      z2 ==> x2 + x1/5 <= 9
    #      z1 + z2 >= 1

    MOI.empty!(model)
    @test MOI.is_empty(model)

    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, MOI.LessThan{Float64}})
    x1 = MOI.add_variable(model)
    x2 = MOI.add_variable(model)
    z1 = MOI.add_variable(model)
    z2 = MOI.add_variable(model)
    MOI.add_constraint(model, z1, MOI.ZeroOne())
    MOI.add_constraint(model, z2, MOI.ZeroOne())
    f1 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0]
    )
    iset1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(8.0))
    MOI.add_constraint(model, f1, iset1)

    f2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(0.2, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0],
    )
    iset2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(9.0))

    MOI.add_constraint(model, f2, iset2)

    # Additional regular constraint.
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x1), MOI.ScalarAffineTerm(1.0, x2)], 0.0),
        MOI.LessThan(10.0),
    )

    # Disjunction z1 ⋁ z2
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, z1), MOI.ScalarAffineTerm(1.0, z2)], 0.0),
        MOI.GreaterThan(1.0),
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0], [x1, x2]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        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.ObjectiveValue()) ≈ 28.75 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x1) ≈ 1.25 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x2) ≈ 8.75 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z1) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z2) ≈ 1.0 atol=atol rtol=rtol
    end
end

function indicator2_test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # linear problem with indicator constraint
    # max  2x1 + 3x2 - 30 z2
    # s.t. x1 + x2 <= 10
    #      z1 ==> x2 <= 8
    #      z2 ==> x2 + x1/5 <= 9
    #      z1 + z2 >= 1

    MOI.empty!(model)
    @test MOI.is_empty(model)

    # This is the same model as indicator_test1, except that the penalty on z2 forces z1 to be 1.

    x1 = MOI.add_variable(model)
    x2 = MOI.add_variable(model)
    z1  = MOI.add_variable(model)
    z2  = MOI.add_variable(model)
    MOI.add_constraint(model, z1, MOI.ZeroOne())
    MOI.add_constraint(model, z2, MOI.ZeroOne())
    f1 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0]
    )
    iset1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(8.0))
    MOI.add_constraint(model, f1, iset1)

    f2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(0.2, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0],
    )
    iset2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(9.0))

    MOI.add_constraint(model, f2, iset2)

    # additional regular constraint
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x1), MOI.ScalarAffineTerm(1.0, x2)], 0.0),
        MOI.LessThan(10.0),
    )

    # disjunction z1 ⋁ z2
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, z1), MOI.ScalarAffineTerm(1.0, z2)], 0.0),
        MOI.GreaterThan(1.0),
    )

    # objective penalized on z2
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0, -30.0], [x1, x2, z2]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        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.ObjectiveValue()) ≈ 28.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x1) ≈ 2.0  atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x2) ≈ 8.0  atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z1) ≈ 1.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z2) ≈ 0.0 atol=atol rtol=rtol
    end
end

function indicator3_test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # linear problem with indicator constraint
    # similar to indicator1_test with reversed z1
    # max  2x1 + 3x2
    # s.t. x1 + x2 <= 10
    #      z1 == 0 ==> x2 <= 8
    #      z2 == 1 ==> x2 + x1/5 <= 9
    #      (1-z1) + z2 >= 1 <=> z2 - z1 >= 0

    MOI.empty!(model)
    @test MOI.is_empty(model)

    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, MOI.LessThan{Float64}})
    x1 = MOI.add_variable(model)
    x2 = MOI.add_variable(model)
    z1 = MOI.add_variable(model)
    z2 = MOI.add_variable(model)
    vc1 = MOI.add_constraint(model, z1, MOI.ZeroOne())
    @test vc1.value == z1.value
    vc2 = MOI.add_constraint(model, z2, MOI.ZeroOne())
    @test vc2.value == z2.value
    f1 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0]
    )
    iset1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.LessThan(8.0))
    MOI.add_constraint(model, f1, iset1)

    f2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(0.2, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0],
    )
    iset2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(9.0))

    MOI.add_constraint(model, f2, iset2)

    # Additional regular constraint.
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x1), MOI.ScalarAffineTerm(1.0, x2)], 0.0),
        MOI.LessThan(10.0),
    )

    # Disjunction (1-z1) ⋁ z2
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, z1), MOI.ScalarAffineTerm(1.0, z2)], 0.0),
        MOI.GreaterThan(0.0),
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0], [x1, x2]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        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.ObjectiveValue()) ≈ 28.75 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x1) ≈ 1.25 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x2) ≈ 8.75 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z1) ≈ 1.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z2) ≈ 1.0 atol=atol rtol=rtol
    end
end

function indicator4_test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # equivalent to indicator1_test with left-hand-side partially in LHS constant
    # linear problem with indicator constraint and
    # max  2x1 + 3x2
    # s.t. x1 + x2 <= 10
    #      z1 ==> x2 - 1 <= 7
    #      z2 ==> x2 + x1/5 + 1 <= 10
    #      z1 + z2 >= 1

    MOI.empty!(model)
    @test MOI.is_empty(model)

    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, MOI.LessThan{Float64}})
    x1 = MOI.add_variable(model)
    x2 = MOI.add_variable(model)
    z1 = MOI.add_variable(model)
    z2 = MOI.add_variable(model)
    MOI.add_constraint(model, z1, MOI.ZeroOne())
    MOI.add_constraint(model, z2, MOI.ZeroOne())
    f1 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, -1.0]
    )
    iset1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(7.0))
    MOI.add_constraint(model, f1, iset1)

    f2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(0.2, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 1.0],
    )
    iset2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(10.0))

    MOI.add_constraint(model, f2, iset2)

    # Additional regular constraint.
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x1), MOI.ScalarAffineTerm(1.0, x2)], 0.0),
        MOI.LessThan(10.0),
    )

    # Disjunction z1 ⋁ z2
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, z1), MOI.ScalarAffineTerm(1.0, z2)], 0.0),
        MOI.GreaterThan(1.0),
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0], [x1, x2]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        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.ObjectiveValue()) ≈ 28.75 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x1) ≈ 1.25 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x2) ≈ 8.75 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z1) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), z2) ≈ 1.0 atol=atol rtol=rtol
    end
end

function _semitest(model::MOI.ModelLike, config::TestConfig{T}, int::Bool) where T
    atol = config.atol
    rtol = config.rtol

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    if !int
        @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Semicontinuous{T})
    else
        @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.Semiinteger{T})
    end

    # 2 variables
    #
    # min  x
    # st   x >= y
    #      if !int
    #           x ∈ {0.0} U [2.0,3.0]
    #      if int
    #           x ∈ {0.0} U {2.0} U {3.0}
    #      y = 0.0

    MOI.empty!(model)
    @test MOI.is_empty(model)

    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    if !int
        vc1 = MOI.add_constraint(model, MOI.SingleVariable(v[1]), MOI.Semicontinuous(T(2), T(3)))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Semicontinuous{T}}()) == 1
    else
        vc1 = MOI.add_constraint(model, MOI.SingleVariable(v[1]), MOI.Semiinteger(T(2), T(3)))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Semiinteger{T}}()) == 1
    end

    vc2 = MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.EqualTo(zero(T)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.EqualTo{T}}()) == 1

    cf = MOI.ScalarAffineFunction{T}(MOI.ScalarAffineTerm{T}.([one(T), -one(T)], v), zero(T))
    c = MOI.add_constraint(model, cf, MOI.GreaterThan(zero(T)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}}()) == 1

    objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.0], v), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.PrimalStatus()) in [ MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT ]
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0, 0] atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ObjectiveBound()) <= 0.0 + atol
    end

    # Change y fixed value

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(one(T)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.0, 1.0] atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0 + atol
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(2)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.0, 2.0] atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0 + atol
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(5//2)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        if !int
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.5 atol=atol rtol=rtol
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.5, 2.5] atol=atol rtol=rtol
            @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol
            @test MOI.get(model, MOI.ObjectiveBound()) <= 2.5 + atol
        else
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [3.0, 2.5] atol=atol rtol=rtol
            @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.5 atol=atol rtol=rtol
            @test MOI.get(model, MOI.ObjectiveBound()) <= 3.0 + atol
        end
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(3)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [3.0, 3.0] atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ObjectiveBound()) <= 3.0 + atol
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(4)))

     if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE ||
            MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE_OR_UNBOUNDED
    end
end

semiconttest(model::MOI.ModelLike, config::TestConfig) = _semitest(model, config, false)
semiinttest(model::MOI.ModelLike, config::TestConfig) = _semitest(model, config, true)

const intlineartests = Dict("knapsack" => knapsacktest,
                            "int1"     => int1test,
                            "int2"     => int2test,
                            "int3"     => int3test,
                            "indicator1"  => indicator1_test,
                            "indicator2"  => indicator2_test,
                            "indicator3"  => indicator3_test,
                            "indicator4"  => indicator4_test,
                            "semiconttest" => semiconttest,
                            "semiinttest" => semiinttest
                           )

@moitestset intlinear


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_reserve_margin(path::AbstractString, setup::Dict, EP::Model)
	temp_ResMar = dual.(EP[:cCapacityResMargin])
	if setup["ParameterScale"] == 1
		temp_ResMar = temp_ResMar * ModelScalingFactor # Convert from MillionUS$/GWh to US$/MWh
	end
	dfResMar = DataFrame(temp_ResMar, :auto)
	CSV.write(joinpath(path, "ReserveMargin.csv"), dfResMar)
	return dfResMar
end

ext: jl

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using Images, MXNet, MLDatasets, StatsBase

DATASET_DIR = "data/101_ObjectCategories"
CLASSES = readdir(DATASET_DIR)

images = Tuple{String,Int}[]

for dir_id = 1:length(CLASSES)

    dir = CLASSES[dir_id]
    all_files = readdir(joinpath(DATASET_DIR, dir))
    keep_files = filter(x -> endswith(x, ".jpg"), all_files)

    append!(images, map(x -> (joinpath(DATASET_DIR, dir, x), dir_id), keep_files))
end

### LOAD THE NEURAL NETWORK
nnet = mx.load_checkpoint("weights/mobilenet-v2/MobiletNet-FE", 0, mx.FeedForward; context = mx.gpu());
nnet.arch = @mx.chain nnet.arch => mx.FullyConnected(num_hidden=length(CLASSES)) => mx.SoftmaxOutput(mx.Variable(:label))

### SHUFFLE DATASET
images = images[shuffle(1:length(images))];

### DEFINE SIZE OF TRAIN AND TEST DATASETS
batch_size = 1000
train_indices = 1:8000
valid_indices = 8000:length(images)

optimizer = mx.ADAM()

for i = 1:100

    println(i, " - ", i+batch_size-1);

    train_indices_batch = sample(train_indices, batch_size, replace = false);
    valid_indices_batch = sample(valid_indices, batch_size, replace = false);

    mx_data_train = mx.zeros((224, 224, 3, batch_size));
    mx_data_train_y = zeros(Int, batch_size);

    mx_data_valid = mx.zeros((224, 224, 3, batch_size));
    mx_data_valid_y = zeros(Int, batch_size);

    for idx = 1:batch_size

        mx_data_train[idx:idx] = reshape(Float16.(channelview(RGB.(imresize(load(images[train_indices_batch[idx]][1]), (224, 224))))), (224, 224, 3, 1));
        mx_data_train_y[idx] = images[train_indices_batch[idx]][2];

        mx_data_valid[idx:idx] = reshape(Float16.(channelview(RGB.(imresize(load(images[valid_indices_batch[idx]][1]), (224, 224))))), (224, 224, 3, 1));
        mx_data_valid_y[idx] = images[valid_indices_batch[idx]][2];
    end

    mx_data_train *= 256.0;
    mx_data_train -= 128.;
    mx_data_train /= 128.;

    mx_data_valid *= 256.0;
    mx_data_valid -= 128.;
    mx_data_valid /= 128.;

    train_data_provider = mx.ArrayDataProvider(:data => mx_data_train, :label => mx_data_train_y, batch_size = 25);
    validation_data_provider = mx.ArrayDataProvider(:data => mx_data_valid, :label => mx_data_valid_y, batch_size = 25);

    mx.fit(nnet, optimizer, train_data_provider, eval_data = validation_data_provider, n_epoch = 5, callbacks = [mx.speedometer()]);

end

### PREPARING THE INPUT
batch_size = 100
train_indices = 1:40000
validation_indices = 40001:50000
optimizer = mx.ADAM();

for i = 1:10

    println(i, " - ", i+batch_size-1)

    train_indices_batch = sample(train_indices, batch_size, replace = false)
    validation_indices_batch = sample(validation_indices, batch_size, replace = false)

    mx_data_train = mx.zeros((224, 224, 3, batch_size));
    mx_data_validation = mx.zeros((224, 224, 3, batch_size));

    map(x -> mx_data_train[x:x] = reshape(Float16.(channelview(imresize(train_x[:, :, :, train_indices_batch[x]], (224, 224)))), (224, 224, 3, 1)), 1:batch_size);
    map(x -> mx_data_validation[x:x] = reshape(Float16.(channelview(imresize(train_x[:, :, :, validation_indices_batch[x]], (224, 224)))), (224, 224, 3, 1)), 1:batch_size);

    mx_data_train *= 256.0;
    mx_data_train -= 128.;
    mx_data_train /= 128.;

    mx_data_validation *= 256.0;
    mx_data_validation -= 128.;
    mx_data_validation /= 128.;

    train_data_provider = mx.ArrayDataProvider(:data => mx_data_train, :label => train_y[train_indices_batch], batch_size = 50);
    validation_data_provider = mx.ArrayDataProvider(:data => mx_data_validation, :label => train_y[validation_indices_batch], batch_size = 100);

    mx.fit(nnet, optimizer, train_data_provider, eval_data = validation_data_provider, n_epoch = 5, callbacks = [mx.speedometer()]);
end

features_train = features[:, 1:4000]
features_validate = features[:, 4001:end];
features_train_y = outputs[1:4000]
features_validate_y = outputs[4001:end]

train_length = 4000
validation_length = 1000

train_data_array  = mx.zeros((size(features_train, 1)..., train_length...));
train_label_array = mx.zeros(train_length);

validation_data_array  = mx.zeros((size(features_train, 1)..., validation_length...));
validation_label_array = mx.zeros(validation_length);

# The number of records we send to the training should be at least number of outcome
for idx = 1:train_length
    train_data_array[idx:idx] = reshape(features_train[:, idx], (size(features_train, 1)..., 1...))
    train_label_array[idx:idx] = features_train_y[idx]
end

for idx = 1:validation_length
    validation_data_array[idx:idx] = reshape(features_validate[:, idx], (size(features_validate, 1)..., 1...))
    validation_label_array[idx:idx] = features_validate_y[idx]
end

train_data_provider = mx.ArrayDataProvider(:data => features_train, :label => train_label_array, batch_size = 500);
validation_data_provider = mx.ArrayDataProvider(:data => features_validate, :label => validation_label_array, batch_size = 500);

arch = @mx.chain mx.Variable(:data) =>
  mx.FullyConnected(num_hidden=128) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=64) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=10) =>
  mx.SoftmaxOutput(mx.Variable(:label))

arch = @mx.chain mx.Variable(:data) =>
  mx.FullyConnected(num_hidden=128) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=64) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=10) =>
  mx.Activation(act_type=:relu) =>
  mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.gpu())

mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = validation_data_provider, n_epoch = 250, callbacks = [mx.speedometer()]);
mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = train_data_provider, n_epoch = 250, callbacks = [mx.speedometer()]);


text:
"""
Solve system of linear equations: L*x = b, where L is Laplacian matrix
"""
function linsolve_cg_v2!( LF::LF3dGrid, x::Array{Float64,1}, b::Array{Float64,1};
                          NiterMax=1000, verbose=false, TOL=5.e-10 )
  #
  Npoints = size(x)[1]
  #
  L_x = apply_Laplacian( LF, x )
  #
  r = b[:] - L_x[:]
  z = r[:]
  p = z[:]

  r_old = r[:]
  z_old = z[:]
  p_old = p[:]
  for iter = 1 : NiterMax
    #
    L_x = apply_Laplacian( LF, p )
    #
    α = dot( r, z ) / dot( p, L_x )
    #
    x[:] = x[:] + α * p[:]  # FIXME use x[:] to force x to be copied, not referenced
    r[:] = r[:] - α * L_x[:]
    #
    norm_res = sqrt( dot( r, r ) )
    if verbose
      @printf("%8d %20.10f\n", iter, norm_res)
    end
    #
    if norm_res < TOL
      if verbose
        @printf("#CONV in linsolve_cg_v2: iter, norm_res: %8d %10.5e\n", iter, norm_res)
      end
      break
    end
    #
    z = r[:]  #XXX call appropriate preconditioner here
    # Polak-Ribiere
    β = dot( z, r - r_old ) / dot( z_old, r_old )
    p = z[:] + β * p_old[:]
    #
    z_old = z[:]
    r_old = r[:]
    p_old = p[:]
  end
  #
  return
  #
end # of function


text:
include("HelperFunctions.jl")

"""
    This function applies a Kalman filter and fixed-interval smoother. The script uses the following model:
        y_t = C_t * Z_t + e_t, for e_t ~ N(0, R)
        Z_t = A * Z_{t-1} + mu_t, for mu_t ~ N(0, Q)
    It then applies a fixed-interval smoother using the results from the KF.
    Throughout this file,
        'm' denotes the number of elements in the state vector Z_t.
        'k' denotes the number of elements (observed variables) in y_t.
        'nobs' denotes the number of time periods for which data are observed.

     parameters:
        y_est : Array
            k x nobs matrix of input data
        A : Array
            m x m transition matrix
        C: Array
            k x m measurement matrix
        Q: Array
            m x m covariance matrix for transition equation residuals (mu_t)
        R: Array
            k x k covariance for measurement matrix residuals (e_t)
        Z0: Array
            1 x m vector, initial value of state
        V0 : Array
            m x m matrix, initial value of factor covariance matrix
    returns: Dict
       Zsmooth => k-by-(nobs+1) matrix, smoothed factor estimates (i.e. Zsmooth[, t + 1] = Z_t|T)
       Vsmooth => k-by-k-by-(nobs+1) array, smoothed factor covariance matrices (i.e. Vsmooth[, , t + 1) = Cov(Z_t|T))
       VmU => Filtered factor posterior covariance
       VVsmooth => k-by-k-by-nobs array, lag 1 factor covariance matrices (i.e. Cov(Z_t, Z_t-1|T))
       loglik => scalar, log-likelihood
"""
function kalman_filter(y_est, A, C, Q, R, Z0, V0)
    # initialize
    m = size(C)[2]
    n_obs = size(y_est)[2]
    Zm = Array{Union{Missing, Float64},2}(undef, m, n_obs) # Z_t | t-1 (prior)
    Vm = Array{Union{Missing, Float64},3}(undef, m, m, n_obs) # V_t | t-1 (prior)
    ZmU = Array{Union{Missing, Float64},2}(undef, m, n_obs+1) # Z_t | t (posterior/updated)
    VmU = Array{Union{Missing, Float64},3}(undef, m, m, n_obs+1) # V_t | t (posterior/updated)
    ZmT = [0.0 for i in 1:(m * (n_obs + 1))] |> x-> reshape(x, m, n_obs + 1) # Z_t | T (smoothed states)
    VmT = [0.0 for i in 1:(m * m * (n_obs + 1))] |> x-> reshape(x, m, m, n_obs + 1) # V_t | T = Cov(Z_t|T) (smoothed factor covariance)
    VmT_lag = Array{Union{Missing, Float64},3}(undef, m, m, n_obs) # Cov(Z_t, Z_t-1|T) (smoothed lag 1 factor covariance)
    loglik = 0

    # initial values
    Zu = Z0 # Z_0|0 (In loop, Zu gives Z_t | t)
    Vu = V0 # V_0|0 (In loop, Vu gives V_t | t)
    ZmU[:, 1] = Zu
    VmU[:,:, 1] = Vu

    # kalman filter
    for t in 1:n_obs
        ### Calculate prior distribution
        Z = A * Zu # Use transition equation to create prior estimate for factor, i.e. Z = Z_t|t-1
        V = A * Vu * transpose(A) + Q # Prior covariance matrix of Z (i.e. V = V_t|t-1)
        V = 0.5 * (V + transpose(V))  # Trick to make symmetric

        ### Calculate posterior distribution
        non_missing = findall(!ismissing, y_est[:,t]) # Remove missing series: These are removed from Y, C, and R
        global yt = y_est[:,t][non_missing]
        global Ct = C[non_missing,:]
        Rt = R[non_missing,non_missing]

        # Check if yt contains no data; if this is the case, replace Zu and Vu with prior
        if length(yt) == 0
            Zu = Z
            Vu = V
        else
            # Steps for variance and population regression coefficients:
            VC = V * transpose(Ct)
            iF = inv(Ct * VC + Rt)
            global VCF = VC * iF # Matrix of population regression coefficients (QuantEcon eqn #4)
            innov = yt - Ct * Z # Gives difference between actual and predicted measurement matrix values
            Zu  = Z  + VCF * innov # Update estimate of factor values (posterior)
            Vu = V - VCF * transpose(VC) # Update covariance matrix (posterior) for time t
            Vu = 0.5 * (Vu + transpose(Vu)) # Trick to make symmetric
            loglik = loglik + 0.5 * (log(det(iF)) - transpose(innov) * iF * innov) # Update log likelihood
        end

        ### Store output
        # Store covariance and observation values for t-1 (priors)
        Zm[:,t] = Z
        Vm[:,:,t] = V
        # Store covariance and state values for t (posteriors), i.e. Zu = Z_t|t & Vu = V_t|t
        ZmU[:,t + 1] = Zu
        VmU[:,:,t + 1] = Vu
    end

    # Store Kalman gain k_t
    if length(yt) == 0
        k_t = zeros(m, m)
    else
        k_t = VCF * Ct
    end

    ### Apply fixed interval smoother
    # Fill the final period of ZmT & VmT with posterior values from KF
    ZmT[:,n_obs + 1] = ZmU[:,n_obs + 1]
    VmT[:,:,n_obs + 1] = VmU[:,:,n_obs + 1]
    VmT_lag[:,:,n_obs] = (I(m) - k_t) * A * VmU[:,:,n_obs] # Initialize VmT_1 lag 1 covariance matrix for final period
    J_2 = VmU[:,:, n_obs] * transpose(A) * pinv(Vm[:,:,n_obs]) # Used for recursion process, see companion file for details

    ### Run smoothing algorithm
    # Loop through time reverse-chronologically (starting at final period nobs)
    for t in n_obs:-1:1
        # Store posterior and prior factor covariance values
        VmUt = VmU[:,:,t]
        Vmt = Vm[:,:,t]
        # Store previous period smoothed factor covariance and lag-1 covariance
        Vt = VmT[:,:,t+1]
        Vt_lag = VmT_lag[:,:,t]
        J_1 = copy(J_2)
        ZmT[:,t] = ZmU[:,t] + J_1 * (ZmT[:,t+1] - A * ZmU[:,t])  # Update smoothed factor estimate
        VmT[:,:,t] = VmUt + J_1 * (Vt - Vmt) * transpose(J_1) # Update smoothed factor covariance matrix
        if t > 1
            J_2 = VmU[:,:,t-1] * transpose(A) * pinv(Vm[:,:,t-1]) # Update weight
            VmT_lag[:,:,t-1] = VmUt * transpose(J_2) + J_1 * (Vt_lag - A * VmUt) * transpose(J_2) # Update lag 1 factor covariance matrix
        end
    end

    ### output
    return Dict(
        :Zsmooth => ZmT,
        :Vsmooth => VmT,
        :VmU => VmU,
        :VVsmooth => VmT_lag,
        :loglik => loglik
    )
end

"""
    This function applies a Kalman filter for news calculation step, when model parameters are already estimated. This procedure only smoothes and fills missing data for a given data matrix
     parameters:
        data : Array
            input data matrix
        output_dfm : Dict
            output of estimate_dfm function
        lag: Int
            number of lags
    returns: Dict
        Plag => Smoothed factor covariance for transition matrix
        Vsmooth => Smoothed factor covariance matrix
        X_smooth => Smoothed data matrix
        F => Smoothed factors
"""
function kalman_filter_constparams(data; output_dfm, lag)
    # initialization
    Z0 = output_dfm[:Z0]; V0 = output_dfm[:V0]; A = output_dfm[:A]; C = output_dfm[:C]; Q = output_dfm[:Q]; R = output_dfm[:R]; means = output_dfm[:means]; stds = output_dfm[:stds];
    stds = output_dfm[:stds]
    means = output_dfm[:means]
    y = transpose(Array(standardize(data)))
    # y = standardize(data) |> j-> j[[sum(.!ismissing.(Array(x))) > 0 for x in eachrow(j)], :] |> Array |> transpose # temporary, should I removing missing before here?
    kalman_output = kalman_filter(y, A, C, Q, R, Z0, V0)

    Vs = kalman_output[:Vsmooth][:,:,1:end-1] # Smoothed factor covariance for transition matrix
    Vf = kalman_output[:VmU][:,:,1:end-1] # Filtered factor posterior covariance
    Zsmooth = kalman_output[:Zsmooth] # Smoothed factors
    Vsmooth = kalman_output[:Vsmooth] # Smoothed covariance value
    Plag = []
    push!(Plag, Vs)

    if lag > 0
        for jk in 1:lag
            tmp_plag = Array{Union{Missing, Float64},3}(undef, size(C)[2], size(C)[2], size(y)[2])
            for jt in size(Plag[1])[3]:-1:(lag + 1)
                As = Vf[:,:, jt - jk] * transpose(A) * pinv(A * Vf[:,:,jt-jk] * transpose(A) + Q)
                tmp_plag[:,:,jt] = As * Plag[jk][:,:,jt]
            end
            push!(Plag, tmp_plag)
        end
    end
    Zsmooth = transpose(Zsmooth)
    x_sm = Zsmooth[1:end-1,:] * transpose(C) # Factors to series representation
    X_smooth = repeat(stds, size(y)[2]) .* x_sm .+ repeat(means, size(y)[2]) # standardized to unstandardized

    return Dict(
        :Plag => Plag,
        :Vsmooth => Vsmooth,
        :X_smooth => X_smooth,
        :F => Zsmooth[1:end-1,:]
    )
end


text:
"""
    Norm

Abstract type for all norms.
"""
abstract type Norm end

# fallback method for operators

LinearAlgebra.opnorm(A::LinearOperator, d::Norm) = opnorm(A, d, d)

function LinearAlgebra.opnorm(A::LinearOperator, d_domain::Norm, d_codomain::Norm)
    codomain_A = codomain(A)
    A_ = coefficients(A)
    @inbounds v₁ = norm(Sequence(codomain_A, view(A_, :, 1)), d_codomain)
    len = size(A, 2)
    v = Matrix{typeof(v₁)}(undef, 1, len)
    @inbounds v[1] = v₁
    @inbounds for i ∈ 2:len
        v[i] = norm(Sequence(codomain_A, view(A_, :, i)), d_codomain)
    end
    return opnorm(LinearOperator(domain(A), ParameterSpace(), v), d_domain)
end

# Sequence spaces

# decay rates

abstract type Weights end

struct GeometricWeights{T<:Real} <: Weights
    rate :: T
    function GeometricWeights{T}(rate::T) where {T<:Real}
        isfinite(rate) & (rate ≥ 1) || return throw(DomainError)
        return new{T}(rate)
    end
end
GeometricWeights(rate::T) where {T<:Real} = GeometricWeights{T}(rate)
rate(weights::GeometricWeights) = weights.rate
weight(weights::GeometricWeights, i::Int) = weights.rate ^ abs(i)
weight(weights::GeometricWeights{<:Interval}, i::Int) = pow(weights.rate, abs(i))

struct AlgebraicWeights{T<:Real} <: Weights
    rate :: T
    function AlgebraicWeights{T}(rate::T) where {T<:Real}
        isfinite(rate) & (rate ≥ 0) || return throw(DomainError)
        return new{T}(rate)
    end
end
AlgebraicWeights(rate::T) where {T<:Real} = AlgebraicWeights{T}(rate)
rate(weights::AlgebraicWeights) = weights.rate
weight(weights::AlgebraicWeights, i::Int) = (1 + abs(i)) ^ weights.rate
weight(weights::AlgebraicWeights{<:Interval}, i::Int) = pow(Interval(1 + abs(i)), weights.rate)

rate(weights::Tuple{Vararg{Weights}}) = map(rate, weights)
@generated function weight(weights::NTuple{N,Weights}, α::NTuple{N,Int}) where {N}
    p = :(weight(weights[1], α[1]))
    for i ∈ 2:N
        p = :(weight(weights[$i], α[$i]) * $p)
    end
    return p
end

function geometricweights(a::Sequence{<:BaseSpace})
    rate = _geometric_rate(space(a), coefficients(a))
    return GeometricWeights(max(one(rate), rate))
end
function geometricweights(a::Sequence{<:TensorSpace})
    rate = _geometric_rate(space(a), coefficients(a))
    return GeometricWeights.(max.(one.(rate), rate))
end

function algebraicweights(a::Sequence{<:BaseSpace})
    rate = _algebraic_rate(space(a), coefficients(a))
    return AlgebraicWeights(max(zero(rate), rate))
end
function algebraicweights(a::Sequence{<:TensorSpace})
    rate = _algebraic_rate(space(a), coefficients(a))
    return AlgebraicWeights.(max.(zero.(rate), rate))
end

_geometric_rate(s::BaseSpace, A::AbstractVector{<:Interval}) = Interval(exp(-_linear_regression(s, log.(mag.(A)))))
_geometric_rate(s::BaseSpace, A::AbstractVector{<:Complex{<:Interval}}) = Interval(exp(-_linear_regression(s, log.(mag.(A)))))
_geometric_rate(s::BaseSpace, A::AbstractVector) = exp(-_linear_regression(s, log.(abs.(A))))
_geometric_rate(s::TensorSpace, A::AbstractVector{<:Interval}) = Interval.(exp.((-).(_linear_regression(s, log.(mag.(A))))))
_geometric_rate(s::TensorSpace, A::AbstractVector{<:Complex{<:Interval}}) = Interval.(exp.((-).(_linear_regression(s, log.(mag.(A))))))
_geometric_rate(s::TensorSpace, A::AbstractVector) = exp.((-).(_linear_regression(s, log.(abs.(A)))))

_algebraic_rate(s::BaseSpace, A::AbstractVector{<:Interval}) = Interval(-_log_linear_regression(s, log.(mag.(A))))
_algebraic_rate(s::BaseSpace, A::AbstractVector{<:Complex{<:Interval}}) = Interval(-_log_linear_regression(s, log.(mag.(A))))
_algebraic_rate(s::BaseSpace, A::AbstractVector) = -_log_linear_regression(s, log.(abs.(A)))
_algebraic_rate(s::TensorSpace, A::AbstractVector{<:Interval}) = Interval.((-).(_log_linear_regression(s, log.(mag.(A)))))
_algebraic_rate(s::TensorSpace, A::AbstractVector{<:Complex{<:Interval}}) = Interval.((-).(_log_linear_regression( s, log.(mag.(A)))))
_algebraic_rate(s::TensorSpace, A::AbstractVector) = (-).(_log_linear_regression(s, log.(abs.(A))))

function _linear_regression(s::TensorSpace{<:NTuple{N,BaseSpace}}, A) where {N}
    A_ = filter(isfinite, A)
    x = ones(Int, length(A_), N+1)
    @inbounds for (i, α) ∈ enumerate(indices(s))
        if isfinite(A[i])
            view(x, i, 2:N+1) .= abs.(α) .+ 1
        end
    end
    x_T = transpose(x)
    r = (x_T * x) \ x_T * A_
    return @inbounds ntuple(i -> r[i+1], Val(N))
end

function _log_linear_regression(s::TensorSpace{<:NTuple{N,BaseSpace}}, A) where {N}
    A_ = filter(isfinite, A)
    x = ones(Float64, length(A_), N+1)
    @inbounds for (i, α) ∈ enumerate(indices(s))
        if isfinite(A[i])
            view(x, i, 2:N+1) .= log.(abs.(α) .+ 1)
        end
    end
    x_T = transpose(x)
    r = (x_T * x) \ x_T * A_
    return @inbounds ntuple(i -> r[i+1], Val(N))
end

# Taylor

function _linear_regression(::Taylor, A)
    n = sum_x = 0
    u = t = sum_A = zero(eltype(A))
    for (i, Aᵢ) ∈ enumerate(A)
        if isfinite(Aᵢ)
            sum_x += i
            u += i*Aᵢ
            t += i*i
            sum_A += Aᵢ
            n += 1
        end
    end
    x̄ = sum_x/n
    r = (u - x̄*sum_A)/(t - sum_x*x̄)
    return ifelse(isfinite(r), r, zero(r))
end

function _log_linear_regression(::Taylor, A)
    n = 0
    sum_x = u = t = sum_A = zero(promote_type(Float64, eltype(A)))
    for (i, Aᵢ) ∈ enumerate(A)
        if isfinite(Aᵢ)
            log_i = log(i)
            sum_x += log_i
            u += log_i*Aᵢ
            t += log_i*log_i
            sum_A += Aᵢ
            n += 1
        end
    end
    x̄ = sum_x/n
    r = (u - x̄*sum_A)/(t - sum_x*x̄)
    return ifelse(isfinite(r), r, zero(r))
end

# Fourier

function _linear_regression(s::Fourier, A)
    ord = order(s)
    n = sum_x = 0
    u = t = sum_A = zero(eltype(A))
    for (i, Aᵢ) ∈ enumerate(A)
        if isfinite(Aᵢ)
            abs_i = abs(i-ord-1)+1
            sum_x += abs_i
            u += abs_i*Aᵢ
            t += abs_i*abs_i
            sum_A += Aᵢ
            n += 1
        end
    end
    x̄ = sum_x/n
    r = (u - x̄*sum_A)/(t - sum_x*x̄)
    return ifelse(isfinite(r), r, zero(r))
end

function _log_linear_regression(s::Fourier, A)
    ord = order(s)
    n = 0
    sum_x = u = t = sum_A = zero(promote_type(Float64, eltype(A)))
    for (i, Aᵢ) ∈ enumerate(A)
        if isfinite(Aᵢ)
            log_abs_i = log(abs(i-ord-1)+1)
            sum_x += log_abs_i
            u += log_abs_i*Aᵢ
            t += log_abs_i*log_abs_i
            sum_A += Aᵢ
            n += 1
        end
    end
    x̄ = sum_x/n
    r = (u - x̄*sum_A)/(t - sum_x*x̄)
    return ifelse(isfinite(r), r, zero(r))
end

# Chebyshev

function _linear_regression(::Chebyshev, A)
    n = sum_x = 0
    u = t = sum_A = zero(eltype(A))
    for (i, Aᵢ) ∈ enumerate(A)
        if isfinite(Aᵢ)
            sum_x += i
            u += i*Aᵢ
            t += i*i
            sum_A += Aᵢ
            n += 1
        end
    end
    x̄ = sum_x/n
    r = (u - x̄*sum_A)/(t - sum_x*x̄)
    return ifelse(isfinite(r), r, zero(r))
end

function _log_linear_regression(::Chebyshev, A)
    n = 0
    sum_x = u = t = sum_A = zero(promote_type(Float64, eltype(A)))
    for (i, Aᵢ) ∈ enumerate(A)
        if isfinite(Aᵢ)
            log_i = log(i)
            sum_x += log_i
            u += log_i*Aᵢ
            t += log_i*log_i
            sum_A += Aᵢ
            n += 1
        end
    end
    x̄ = sum_x/n
    r = (u - x̄*sum_A)/(t - sum_x*x̄)
    return ifelse(isfinite(r), r, zero(r))
end

# ℓᵖ norm

"""
    ℓᵖNorm{T<:Real} <: Norm

Norm of the ``\\ell^p`` space.

Fields:
- `p :: T`
"""
struct ℓᵖNorm{T<:Real} <: Norm
    p :: T
    function ℓᵖNorm{T}(p::T) where {T<:Real}
        p ≥ 1 || return throw(DomainError(p, "ℓᵖNorm is only defined for real numbers greater than 1 and possibly infinite"))
        return new{T}(p)
    end
end

ℓᵖNorm(p::T) where {T<:Real} = ℓᵖNorm{T}(p)

LinearAlgebra.norm(a::Sequence, p::Real=2) = norm(a, ℓᵖNorm(p))
LinearAlgebra.opnorm(A::LinearOperator, p::Real=2) = opnorm(A, ℓᵖNorm(p))

LinearAlgebra.norm(a::Sequence, d::ℓᵖNorm) = norm(coefficients(a), d.p)

LinearAlgebra.opnorm(A::LinearOperator{<:VectorSpace,ParameterSpace}, d::ℓᵖNorm) =
    opnorm(coefficients(A), d.p)

function LinearAlgebra.opnorm(A::LinearOperator, d_domain::ℓᵖNorm, d_codomain::ℓᵖNorm)
    d_domain.p == d_codomain.p && return opnorm(coefficients(A), d_domain.p)
    v = map(Aᵢ -> norm(Aᵢ, d_codomain.p), eachcol(coefficients(A)))
    return opnorm(transpose(v), d_domain.p)
end

for T ∈ (:Interval, :(Complex{<:Interval}))
    @eval begin
        LinearAlgebra.norm(a::Sequence{<:VectorSpace,<:AbstractVector{<:$T}}, d::ℓᵖNorm) = norm(Interval.(mag.(coefficients(a))), d.p)

        LinearAlgebra.opnorm(A::LinearOperator{<:VectorSpace,ParameterSpace,<:AbstractMatrix{<:$T}}, d::ℓᵖNorm) =
            opnorm(Interval.(mag.(coefficients(A))), d.p)

        function LinearAlgebra.opnorm(A::LinearOperator{<:VectorSpace,<:VectorSpace,<:AbstractMatrix{<:$T}}, d_domain::ℓᵖNorm, d_codomain::ℓᵖNorm)
            mA = Interval.(mag.(coefficients(A)))
            d_domain.p == d_codomain.p && return opnorm(mA, d.p)
            v = map(mAᵢ -> Interval(mag(norm(mAᵢ, d_codomain.p))), eachcol(mA))
            return opnorm(transpose(v), d_domain.p)
        end
    end
end

# weighted ℓ¹ norm

"""
    Weightedℓ¹Norm{T<:Union{Weights,Tuple{Vararg{Weights}}}} <: Norm

Norm of the weighted ``\\ell^1`` space.

Fields:
- `weights :: T`
"""
struct Weightedℓ¹Norm{T<:Union{Weights,Tuple{Vararg{Weights}}}} <: Norm
    weights :: T
end

LinearAlgebra.norm(a::Sequence{<:BaseSpace}, d::Weightedℓ¹Norm{<:Weights}) =
    _apply(d, space(a), coefficients(a))

function LinearAlgebra.norm(a::Sequence{TensorSpace{T}}, d::Weightedℓ¹Norm{<:NTuple{N,Weights}}) where {N,T<:NTuple{N,BaseSpace}}
    space_a = space(a)
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    return _apply(d, space_a, A)
end

_apply(d::Weightedℓ¹Norm, space::TensorSpace, A) =
    @inbounds _apply(Weightedℓ¹Norm(d.weights[1]), space[1], _apply(Weightedℓ¹Norm(Base.tail(d.weights)), Base.tail(space), A))

_apply(d::Weightedℓ¹Norm, space::TensorSpace{<:Tuple{BaseSpace}}, A) =
    @inbounds _apply(Weightedℓ¹Norm(d.weights[1]), space[1], A)

LinearAlgebra.opnorm(A::LinearOperator{<:BaseSpace,ParameterSpace}, d::Weightedℓ¹Norm{<:Weights}) =
    _apply_dual(d, domain(A), vec(coefficients(A)))

function LinearAlgebra.opnorm(A::LinearOperator{TensorSpace{T},ParameterSpace}, d::Weightedℓ¹Norm{<:NTuple{N,Weights}}) where {N,T<:NTuple{N,BaseSpace}}
    domain_A = domain(A)
    A_ = _no_alloc_reshape(coefficients(A), dimensions(domain_A))
    return _apply_dual(d, domain_A, A_)
end

_apply_dual(d::Weightedℓ¹Norm, space::TensorSpace, A) =
    @inbounds _apply_dual(Weightedℓ¹Norm(d.weights[1]), space[1], _apply_dual(Weightedℓ¹Norm(Base.tail(d.weights)), Base.tail(space), A))

_apply_dual(d::Weightedℓ¹Norm, space::TensorSpace{<:Tuple{BaseSpace}}, A) =
    @inbounds _apply_dual(Weightedℓ¹Norm(d.weights[1]), space[1], A)

# Taylor

function _apply(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Taylor, A::AbstractVector)
    ν = rate(d.weights)
    ord = order(space)
    @inbounds s = abs(A[ord+1]) * one(ν)
    @inbounds for i ∈ ord-1:-1:0
        s = s * ν + abs(A[i+1])
    end
    return s
end
function _apply(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Taylor, A::AbstractVector)
    @inbounds s = abs(A[1]) * weight(d.weights, 0)
    @inbounds for i ∈ 1:order(space)
        s += abs(A[i+1]) * weight(d.weights, i)
    end
    return s
end

function _apply(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Taylor, A::AbstractArray{T,N}) where {T,N}
    ν = rate(d.weights)
    CoefType = typeof(abs(zero(T))*ν)
    ord = order(space)
    @inbounds Aᵢ = selectdim(A, N, ord+1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ ord-1:-1:0
        s .= s .* ν .+ abs.(selectdim(A, N, i+1))
    end
    return s
end
function _apply(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Taylor, A::AbstractArray{T,N}) where {T,N}
    CoefType = typeof(abs(zero(T))*weight(d.weights, 0))
    @inbounds Aᵢ = selectdim(A, N, 1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:order(space)
        s .+= abs.(selectdim(A, N, i+1)) .* weight(d.weights, i)
    end
    return s
end

function _apply_dual(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Taylor, A::AbstractVector{T}) where {T}
    ν = rate(d.weights)
    ν⁻¹ = abs(one(T))/ν
    ν⁻ⁱ = one(ν⁻¹)
    @inbounds s = abs(A[1]) * ν⁻ⁱ
    @inbounds for i ∈ 1:order(space)
        ν⁻ⁱ *= ν⁻¹
        s = max(s, abs(A[i+1]) * ν⁻ⁱ)
    end
    return s
end
function _apply_dual(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Taylor, A::AbstractVector{T}) where {T}
    @inbounds s = abs(A[1]) / weight(d.weights, 0)
    @inbounds for i ∈ 1:order(space)
        s = max(s, abs(A[i+1]) / weight(d.weights, i))
    end
    return s
end

function _apply_dual(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Taylor, A::AbstractArray{T,N}) where {T,N}
    ν = rate(d.weights)
    ν⁻¹ = abs(one(T))/ν
    ν⁻ⁱ = one(ν⁻¹)
    CoefType = typeof(ν⁻¹)
    @inbounds Aᵢ = selectdim(A, N, 1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:order(space)
        ν⁻ⁱ *= ν⁻¹
        s .= max.(s, abs.(selectdim(A, N, i+1)) .* ν⁻ⁱ)
    end
    return s
end
function _apply_dual(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Taylor, A::AbstractArray{T,N}) where {T,N}
    CoefType = typeof(abs(zero(T))/weight(d.weights, 0))
    @inbounds Aᵢ = selectdim(A, N, 1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:order(space)
        s .= max.(s, abs.(selectdim(A, N, i+1)) ./ weight(d.weights, i))
    end
    return s
end

# Fourier

function _apply(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Fourier, A::AbstractVector)
    ν = rate(d.weights)
    ord = order(space)
    if ord == 0
        return abs(A[1]) * one(ν)
    else
        @inbounds s = (abs(A[1]) + abs(A[2ord+1])) * one(ν)
        @inbounds for i ∈ ord-1:-1:1
            s = s * ν + abs(A[ord+1-i]) + abs(A[ord+1+i])
        end
        @inbounds s = s * ν + abs(A[ord+1])
        return s
    end
end
function _apply(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Fourier, A::AbstractVector)
    ord = order(space)
    @inbounds s = abs(A[ord+1]) * weight(d.weights, 0)
    @inbounds for i ∈ 1:ord
        s += (abs(A[ord+1-i]) + abs(A[ord+1+i])) * weight(d.weights, i)
    end
    return s
end

function _apply(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Fourier, A::AbstractArray{T,N}) where {T,N}
    ν = rate(d.weights)
    CoefType = typeof(abs(zero(T))*ν)
    ord = order(space)
    @inbounds A₋ᵢ = selectdim(A, N, 1)
    s = Array{CoefType,N-1}(undef, size(A₋ᵢ))
    if ord == 0
        s .= abs.(A₋ᵢ)
    else
        @inbounds s .= abs.(selectdim(A, N, 2ord+1)) .+ abs.(A₋ᵢ)
        @inbounds for i ∈ ord-1:-1:1
            s .= s .* ν .+ abs.(selectdim(A, N, ord+1-i)) .+ abs.(selectdim(A, N, ord+1+i))
        end
        @inbounds s .= s .* ν .+ abs.(selectdim(A, N, ord+1))
    end
    return s
end
function _apply(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Fourier, A::AbstractArray{T,N}) where {T,N}
    CoefType = typeof(abs(zero(T))*weight(d.weights, 0))
    ord = order(space)
    @inbounds Aᵢ = selectdim(A, N, ord+1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    @inbounds s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:ord
        s .+= (abs.(selectdim(A, N, ord+1-i)) .+ abs.(selectdim(A, N, ord+1+i))) .* weight(d.weights, i)
    end
    return s
end

function _apply_dual(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Fourier, A::AbstractVector{T}) where {T}
    ν = rate(d.weights)
    ν⁻¹ = abs(one(T))/ν
    ν⁻ⁱ = one(ν⁻¹)
    ord = order(space)
    @inbounds s = abs(A[ord+1]) * ν⁻ⁱ
    @inbounds for i ∈ 1:ord
        ν⁻ⁱ *= ν⁻¹
        s = max(s, abs(A[ord+1+i]) * ν⁻ⁱ, abs(A[ord+1-i]) * ν⁻ⁱ)
    end
    return s
end
function _apply_dual(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Fourier, A::AbstractVector{T}) where {T}
    ord = order(space)
    @inbounds s = abs(A[ord+1]) / weight(d.weights, 0)
    @inbounds for i ∈ 1:ord
        x = inv(weight(d.weights, i))
        s = max(s, abs(A[ord+1+i]) * x, abs(A[ord+1-i]) * x)
    end
    return s
end

function _apply_dual(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Fourier, A::AbstractArray{T,N}) where {T,N}
    ν = rate(d.weights)
    ν⁻¹ = abs(one(T))/ν
    ν⁻ⁱ = one(ν⁻¹)
    CoefType = typeof(ν⁻¹)
    ord = order(space)
    @inbounds Aᵢ = selectdim(A, N, ord+1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:ord
        ν⁻ⁱ *= ν⁻¹
        s .= max.(s, abs.(selectdim(A, N, ord+1-i)) .* ν⁻ⁱ, abs.(selectdim(A, N, ord+1+i)) .* ν⁻ⁱ)
    end
    return s
end
function _apply_dual(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Fourier, A::AbstractArray{T,N}) where {T,N}
    CoefType = typeof(abs(zero(T))/weight(d.weights, 0))
    ord = order(space)
    @inbounds Aᵢ = selectdim(A, N, ord+1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:ord
        x = inv(weight(d.weights, i))
        s .= max.(s, abs.(selectdim(A, N, ord+1-i)) ./ x, abs.(selectdim(A, N, ord+1+i)) ./ x)
    end
    return s
end

# Chebyshev

function _apply(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Chebyshev, A::AbstractVector)
    ν = rate(d.weights)
    ord = order(space)
    @inbounds s = abs(A[ord+1]) * one(ν)
    @inbounds for i ∈ ord-1:-1:0
        s = s * ν + abs(A[i+1])
    end
    return s
end
function _apply(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Chebyshev, A::AbstractVector)
    @inbounds s = abs(A[1]) * weight(d.weights, 0)
    @inbounds for i ∈ 1:order(space)
        s += abs(A[i+1]) * weight(d.weights, i)
    end
    return s
end

function _apply(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Chebyshev, A::AbstractArray{T,N}) where {T,N}
    ν = rate(d.weights)
    CoefType = typeof(abs(zero(T))*ν)
    ord = order(space)
    @inbounds Aᵢ = selectdim(A, N, ord+1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ ord-1:-1:0
        s .= s .* ν .+ abs.(selectdim(A, N, i+1))
    end
    return s
end
function _apply(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Chebyshev, A::AbstractArray{T,N}) where {T,N}
    CoefType = typeof(abs(zero(T))*weight(d.weights, 0))
    @inbounds Aᵢ = selectdim(A, N, 1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:order(space)
        s .+= abs.(selectdim(A, N, i+1)) .* weight(d.weights, i)
    end
    return s
end

function _apply_dual(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Chebyshev, A::AbstractVector{T}) where {T}
    ν = rate(d.weights)
    ν⁻¹ = abs(one(T))/ν
    ν⁻ⁱ = one(ν⁻¹)
    @inbounds s = abs(A[1]) * ν⁻ⁱ
    @inbounds for i ∈ 1:order(space)
        ν⁻ⁱ *= ν⁻¹
        s = max(s, abs(A[i+1]) * ν⁻ⁱ)
    end
    return s
end
function _apply_dual(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Chebyshev, A::AbstractVector{T}) where {T}
    @inbounds s = abs(A[1]) / weight(d.weights, 0)
    @inbounds for i ∈ 1:order(space)
        s = max(s, abs(A[i+1]) / weight(d.weights, i))
    end
    return s
end

function _apply_dual(d::Weightedℓ¹Norm{<:GeometricWeights}, space::Chebyshev, A::AbstractArray{T,N}) where {T,N}
    ν = rate(d.weights)
    ν⁻¹ = abs(one(T))/ν
    ν⁻ⁱ = one(ν⁻¹)
    CoefType = typeof(ν⁻¹)
    @inbounds Aᵢ = selectdim(A, N, 1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:order(space)
        ν⁻ⁱ *= ν⁻¹
        s .= max.(s, abs.(selectdim(A, N, i+1)) .* ν⁻ⁱ)
    end
    return s
end
function _apply_dual(d::Weightedℓ¹Norm{<:AlgebraicWeights}, space::Chebyshev, A::AbstractArray{T,N}) where {T,N}
    CoefType = typeof(abs(zero(T))/weight(d.weights, 0))
    @inbounds Aᵢ = selectdim(A, N, 1)
    s = Array{CoefType,N-1}(undef, size(Aᵢ))
    s .= abs.(Aᵢ)
    @inbounds for i ∈ 1:order(space)
        s .= max.(s, abs.(selectdim(A, N, i+1)) ./ weight(d.weights, i))
    end
    return s
end

# Hˢ Sobolev norm

"""
    𝐻ˢNorm{T<:Real} <: Norm

Norm of the ``H^s`` Sobolev space.

Fields:
- `exponent :: T`
"""
struct 𝐻ˢNorm{T<:Real} <: Norm
    exponent :: T
    function 𝐻ˢNorm{T}(exponent::T) where {T<:Real}
        isfinite(exponent) & (exponent > 0) || return throw(DomainError(exponent, "𝐻ˢNorm is only defined for real numbers greater than 1"))
        return new{T}(exponent)
    end
end

𝐻ˢNorm(exponent::T) where {T<:Real} = 𝐻ˢNorm{T}(exponent)

LinearAlgebra.norm(a::Sequence{<:BaseSpace}, d::𝐻ˢNorm) =
    _apply(d, space(a), coefficients(a))

function LinearAlgebra.norm(a::Sequence{<:TensorSpace}, d::𝐻ˢNorm)
    space_a = space(a)
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    return _apply(d, space_a, A)
end

_apply(d::𝐻ˢNorm, space::TensorSpace, A) =
    @inbounds _apply(d, space[1], _apply(d, Base.tail(space), A))

_apply(d::𝐻ˢNorm, space::TensorSpace{<:Tuple{BaseSpace}}, A) =
    @inbounds _apply(d, space[1], A)

LinearAlgebra.opnorm(A::LinearOperator{<:BaseSpace,ParameterSpace}, d::𝐻ˢNorm) =
    _apply_dual(d, domain(A), vec(coefficients(A)))

function LinearAlgebra.opnorm(A::LinearOperator{<:TensorSpace,ParameterSpace}, d::𝐻ˢNorm)
    domain_A = domain(A)
    A_ = _no_alloc_reshape(coefficients(A), dimensions(domain_A))
    return _apply_dual(d, domain_A, A_)
end

_apply_dual(d::𝐻ˢNorm, space::TensorSpace, A) =
    @inbounds _apply_dual(d, space[1], _apply_dual(d, Base.tail(space), A))

_apply_dual(d::𝐻ˢNorm, space::TensorSpace{<:Tuple{BaseSpace}}, A) =
    @inbounds _apply_dual(d, space[1], A)

# Fourier

function _apply(d::𝐻ˢNorm, space::Fourier, A::AbstractVector)
    s = d.exponent
    un = one(eltype(A))
    ord = order(space)
    @inbounds x = abs2(A[ord+1]) * (un + 0)^s
    @inbounds for i ∈ 1:ord
        x += (abs2(A[ord+1-i]) + abs2(A[ord+1+i])) * (un + i*i)^s
    end
    return sqrt(x)
end

function _apply(d::𝐻ˢNorm, space::Fourier, A::AbstractArray{T,N}) where {T,N}
    s = d.exponent
    un = one(T)
    CoefType = typeof(sqrt(abs2(zero(T))*(un+0)^s))
    ord = order(space)
    @inbounds Aᵢ = selectdim(A, N, ord+1)
    x = Array{CoefType,N-1}(undef, size(Aᵢ))
    x .= abs.(Aᵢ)
    @inbounds for i ∈ 1:ord
        x .+= (abs2.(selectdim(A, N, ord+1+i)) .+ abs2.(selectdim(A, N, ord+1-i))) .* (un + i*i)^s
    end
    x .= sqrt.(x)
    return x
end

function _apply_dual(d::𝐻ˢNorm, space::Fourier, A::AbstractVector)
    s = d.exponent
    un = one(eltype(A))
    ord = order(space)
    @inbounds x = abs2(A[ord+1]) / (un + 0)^s
    @inbounds for i ∈ 1:ord
        x += (abs2(A[ord+1-i]) + abs2(A[ord+1+i])) / (un + i*i)^s
    end
    return sqrt(x)
end

function _apply_dual(d::𝐻ˢNorm, space::Fourier, A::AbstractArray{T,N}) where {T,N}
    s = d.exponent
    un = one(T)
    CoefType = typeof(sqrt(abs2(zero(T))/(un+0)^s))
    ord = order(space)
    @inbounds Aᵢ = selectdim(A, N, ord+1)
    x = Array{CoefType,N-1}(undef, size(Aᵢ))
    x .= abs2.(Aᵢ)
    @inbounds for i ∈ 1:ord
        x .+= (abs2.(selectdim(A, N, ord+1+i)) .+ abs2.(selectdim(A, N, ord+1-i))) ./ (un + i*i)^s
    end
    x .= sqrt.(x)
    return x
end

# Cartesian spaces

"""
    CartesianPowerNorm{T<:Norm,S<:Norm} <: Norm

Norm of a cartesian space comprised of spaces with the same norm.

Fields:
- `inner :: T`
- `outer :: S`
"""
struct CartesianPowerNorm{T<:Norm,S<:Norm} <: Norm
    inner :: T
    outer :: S
end

function LinearAlgebra.norm(a::Sequence{<:CartesianSpace}, d::CartesianPowerNorm)
    s = CartesianPower(ParameterSpace(), nb_cartesian_product(space(a)))
    v = map(aᵢ -> norm(aᵢ, d.inner), eachcomponent(a))
    return norm(Sequence(s, v), d.outer)
end

function LinearAlgebra.opnorm(A::LinearOperator{<:CartesianSpace,ParameterSpace}, d::CartesianPowerNorm)
    s = CartesianPower(ParameterSpace(), nb_cartesian_product(domain(A)))
    v = map(Aᵢ -> opnorm(Aᵢ, d.inner), eachcomponent(A))
    return opnorm(LinearOperator(s, ParameterSpace(), v), d.outer)
end

"""
    CartesianProductNorm{T<:Tuple{Vararg{Norm}}, S<:Norm} <: Norm

Norm of a cartesian space comprised of spaces with different norms.

Fields:
- `inner :: T`
- `outer :: S`
"""
struct CartesianProductNorm{T<:Tuple{Vararg{Norm}}, S<:Norm} <: Norm
    inner :: T
    outer :: S
end

function LinearAlgebra.norm(a::Sequence{<:CartesianSpace}, d::CartesianProductNorm{<:NTuple{N,Norm}}) where {N}
    n = nb_cartesian_product(space(a))
    n == N || return throw(DimensionMismatch)
    s = CartesianPower(ParameterSpace(), n)
    v = map((aᵢ, dᵢ) -> norm(aᵢ, dᵢ), eachcomponent(a), d.inner)
    return norm(Sequence(s, v), d.outer)
end

function LinearAlgebra.opnorm(A::LinearOperator{<:CartesianSpace,ParameterSpace}, d::CartesianProductNorm{<:NTuple{N,Norm}}) where {N}
    n = nb_cartesian_product(domain(A))
    n == N || return throw(DimensionMismatch)
    s = CartesianPower(ParameterSpace(), n)
    v = map((Aᵢ, dᵢ) -> opnorm(Aᵢ, dᵢ), eachcomponent(A), d.inner)
    return opnorm(LinearOperator(s, ParameterSpace(), transpose(v)), d.outer)
end


text:
using Test
using PyCall
using RobotOS
RobotOS.debug(true)

#Generally, later tests rely on things defined in previous tests, so the order is important
include("rospy.jl")
include("time.jl")
include("typegeneration.jl")
include("pubsub.jl")
include("services.jl")


text:
import Base: length, GMP

export SortPoints, lowerconvexhull, newtonpolygon, Line, Polygon, PartialPolygon

#using Nemo.Fields

mutable struct Line
  points :: Array{Tuple{fmpq, fmpq},1}
  slope :: fmpq

  function Line(points::Array{Tuple{fmpq, fmpq},1})
    line = new()
    line.points = points
    line.slope = slope(points[1],points[2])
    return line
  end

  function Line(a::Tuple{fmpq,fmpq},b::Tuple{fmpq,fmpq})
    return Line([a,b])
  end
end

mutable struct Polygon
  lines :: Array{Line,1}
  slopes :: Array{fmpq ,1}

  function Polygon(lines::Array{Line,1}; sorted = false)
    polygon = new()
    if sorted
      polyon.lines = lines
    else
      for j = 1:length(lines)-1
        iMin = j
        for i = j+1:length(lines)
          if (lines[i]).points[1][1] < lines[iMin].points[1][1]
            iMin = i
          end
        end
        if iMin != j
          lines[j],lines[iMin] = lines[iMin],lines[j]
        end
      end
      polygon.lines = lines
    end

    polygon.slopes = Array{fmpq}(length(lines))

    for i in 1:length(lines)
      polygon.slopes[i] = lines[i].slope
    end
    return polygon
  end
end

function sortpoints(x::Array{Tuple{fmpq,fmpq},1})
  for j = 1:length(x)-1
    iMin = j
    for i = j+1:length(x)
      if x[i][1] < x[iMin][1]
        iMin = i
      end
    end
    if iMin != j
      x[j],x[iMin] = x[iMin],x[j]
    end
  end
  return x
end

function lowerconvexhull(points::Array{Tuple{fmpq, fmpq},1})

  points = sortpoints(points)

  # Take care of trivial case with 1 or 2 elements

  if length(points) == 1
    println("ups")
  elseif length(points) == 2
    return Line(points)
  end

  pointsonconvexhull = [ points[length(points)-1], points[length(points)] ]

  n = length(points)-2

  oldslope = slope(pointsonconvexhull[1],pointsonconvexhull[2])

  for i = 1:n
    newslope = slope(points[n-i+1], pointsonconvexhull[1])
    if newslope >= oldslope
      pointsonconvexhull[1] = points[n-i+1]
    else
      unshift!(pointsonconvexhull,points[n-i+1])
    end
    oldslope = newslope
  end

  t = Line[]
  for i=1:length(pointsonconvexhull)-1
    push!(t,Line(pointsonconvexhull[i],pointsonconvexhull[i+1]))
  end
  return Polygon(t)
end

function slope(a::Tuple{fmpq,fmpq},b::Tuple{fmpq,fmpq})
  return (b[2]-a[2])//(b[1]-a[1])
end

function newtonpolygon(f::Union(fmpz_poly,fmpq_poly),p::fmpz)
  a = Array{Tuple{fmpq, fmpq}}(degree(f)+1)
  if !(isprime(p))
    #throw some error
  end
  for i = 0:degree(f)
    iszero(coeff(f,i)) ? continue : a[i+1] = (i,valuation(coeff(f,i),p))
  end
  return lowerconvexhull(a)
end

# Remove all lines where the slop is <= t

function PartialPolygon(P::Polygon,t::Rational{BigInt})

  function tem(s::Line)
    if s.slope > t
      return true
    end
    return false
  end

  m = findfirst(tem,P.lines)
  println(m)
  return true
end


text:
"""
    array_response(x, y, sx1, sx2, sy1, sy2, ds, f1, f2, df, s_deg=false) -> arf::ArrayResponse
    array_response(x, y, smax, ds, f1, f2, df, s_deg=false) -> arf::ArrayResponse
    array_response(traces::Array{<:Seis.Trace}, args...) -> arf::ArrayResponse
    array_response(stations::Array{<:Seis.Station}, args...) -> arf

Return the array response function `arf.power` (an [`ArrayResponse`](@ref)) at
slownesses with east and north components `arf.sx` and `arf.sy` respectively,
in s/km.  The ARF is computed between frequencies `f1` and `f2` Hz, spaced by `df` Hz.

Provide either explicit limits on slowness `(sx1, sx2)` and `(sy1, sy2)` in the
east and north components respectively, or the maximum slowness `smax`,
plus the slowness spacing `ds`.

`x` and `y` are the station coordinates in km; or an array of `Seis.Trace` `traces`
can be given.

Default input units of slowness are s/km.  Give `s_deg` as `true`
to use s/°.  Note that values in `arf` are always given in s/km.
"""
function array_response(x, y, sx1, sx2, sy1, sy2, ds, f1, f2, df, s_deg=false)
    if s_deg
        sx1, sx2, sy1, sy2, ds = s_per_km.((sx1, sx2, sy1, sy2, ds))
    end
    arf = ArrayResponse{eltype(x)}(x, y, sx1, sx2, sy1, sy2, ds, f1, f2, df)
    _compute_arf!(arf)
end

array_response(x, y, smax, ds, f1, f2, df, s_deg::Bool=false) =
    array_response(x, y, -smax, smax, -smax, smax, ds, f1, f2, df, s_deg)
array_response(s::AbstractArray{<:Seis.Station}, args...) =
    array_response(array_geometry(s)[1:2]..., args...)
array_response(s::TraceArray, args...) =
    array_response(array_geometry(s)[1:2]..., args...)

"""
    _compute_arf!(arf)

Do the actual ARF calculation for an `ArrayResponse`.
"""
function _compute_arf!(arf::ArrayResponse{T}) where T
    @inbounds for (j, slowy) in enumerate(arf.sy), (i, slowx) in enumerate(arf.sx)
        for f in arf.freqs
            ω = T(2)*π*f
            pow = zero(T)*im
            @simd for istat in 1:length(arf.x)
                pow += cis(ω*(slowx*arf.x[istat] + slowy*arf.y[istat]))
            end
            arf.power[i,j] += abs2(pow)
        end
    end
    arf
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
# Solves the Galewsky test case for the shallow water equations on a sphere
# of physical radius 6371220m. Involves a shear flow instability of a zonal
# jet triggered by an initial gravity wave.
# reference:
#   Galewsky, Scott and Polvani (2004) Tellus, 56A 429-440

const H₀ = 10000.0

function uθ(θϕr)
  θ,ϕ,r = θϕr
  ϵ     = 1.0e-8
  umax  = 80.0
  ϕ₁    = π/7
  ϕ₂    = π/2 - ϕ₁
  en    = exp(-4.0/((ϕ₂ - ϕ₁)*(ϕ₂ - ϕ₁)))
  u     = 0.0
  if ϕ > ϕ₁ + ϵ && ϕ < ϕ₂ - ϵ
    u = (umax/en)*exp(1.0/((ϕ - ϕ₁)*(ϕ - ϕ₂)))
  end
  u
end

# Initial velocity
function u₀(xyz)
  θϕr = xyz2θϕr(xyz)
  u   = uθ(θϕr)
  spherical_to_cartesian_matrix(θϕr)⋅VectorValue(u,0,0)
end

# Initial fluid depth
function h₀(xyz)
  θϕr   = xyz2θϕr(xyz)
  x,y,z = xyz
  θ,ϕ,r = θϕr
  h     = H₀
  hh    = 120.0
  α     = 1.0/3.0
  β     = 1.0/15.0
  ϕ₂    = π/4
  ni    = 1000
  ϕₚ    = 0.0
  dϕ    = abs(ϕ/ni)
  sgn   = 1.0
  if ϕ < 0.0
    sgn = -1.0
  end
  for i in 1:ni
    ϕₚ   = ϕₚ + sgn*dϕ
    _θϕr = VectorValue(θ,ϕₚ,r)
    u    = uθ(_θϕr)
    _f   = 2.0*Ωₑ*sin(ϕₚ)
    h    = h - rₑ*u*(_f + tan(ϕₚ)*u/rₑ)*dϕ/g
  end
  h = h + hh*cos(ϕ)*exp(-1.0*(θ/α)*(θ/α))*exp(-1.0*((ϕ₂ - ϕ)/β)*((ϕ₂ - ϕ)/β))
  h
end

# Topography
function topography(xyz)
  0.0
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
@graphs begin
	@node phase_3 = Start() -> p3_strat
	@node p3_strat = CheckStrategy("military") -> [n_true=p3_mili_1, n_false=p3_corr_1]

	# Military
	@node p3_mili_1 = BinaryCondition("The Fellowship is on the Mordor track.") -> [n_true = p3_mili_1_yes, n_false = p3_mili_2]
	@node p3_mili_1_yes = PerformAction("Assign the maximum allowed number of dice to the hunt pool.") -> p3_mili_end_phase

	@node p3_mili_2 = BinaryCondition("The Fellowship's progress is greater than 5.") -> [n_true = p3_mili_2_yes, n_false = p3_mili_3]
	@node p3_mili_2_yes = PerformAction("Assign 2 dice to the hunt pool.") -> p3_mili_end_phase

	@node p3_mili_3 = BinaryCondition("The Fellowship is on the starting position and its progress is 0.") -> [n_true = p3_mili_3_yes, n_false = p3_mili_3_no]
	@node p3_mili_3_yes = PerformAction("Assign 0 dice to the hunt pool.") -> p3_mili_end_phase
	@node p3_mili_3_no = PerformAction("Assign 1 dice to the hunt pool.") -> p3_mili_end_phase

	@node p3_mili_end_phase = End("End of Phase") -> []

	# Corruption
	@node p3_corr_1 = BinaryCondition("The Fellowship is on the starting position and its progess is 0.") -> [n_true = p3_corr_1_yes, n_false = p3_corr_2]
	@node p3_corr_1_yes = PerformAction("Roll a d6. On 4+, assign 1 die to the hunt pool, otherwise do nothing.") -> p3_corr_end_phase

	@node p3_corr_2 = BinaryCondition( " The Fellowship is on the Mordor track.  ") -> [n_true = p3_corr_2_yes, n_false = p3_corr_3]
	@node p3_corr_2_yes = PerformAction( " Assign the maximum allowed number of dice to the hunt pool.  ") -> p3_corr_end_phase

	@node p3_corr_3 = BinaryCondition("""
									  A *mobile* army is adjacent to its *target* which provides enough victory points to win the game.
									  Or, the Shadow have 7 dice.
									  """
									  ) -> [n_true = p3_corr_3_yes, n_false = p3_corr_4]
	@node p3_corr_3_yes = PerformAction( " Assign 1 die to the hunt pool.  ") -> p3_corr_end_phase

	@node p3_corr_4 = BinaryCondition( " The Fellowship's progress is greater than 4.  ") -> [n_true = p3_corr_4_yes, n_false = p3_corr_5]
	@node p3_corr_4_yes = PerformAction( " Assign 2 dice to the hunt pool.  ") -> p3_corr_end_phase

	@node p3_corr_5 = BinaryCondition("""
									  The Fellowship's shortest path to Mordor leads via a Shadow stronghold and the progress allows them to pass it or be within 2 steps from it.
									  """) -> [n_true = p3_corr_5_yes, n_false = p3_corr_5_no]
	@node p3_corr_5_yes = PerformAction("Assign 2 dice to the hunt pool.") -> p3_corr_end_phase
	@node p3_corr_5_no = PerformAction("Assign 1 dice to the hunt pool.") -> p3_corr_end_phase

	@node p3_corr_end_phase = End("End of Phase") -> []
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
# developed with Julia 1.1.1
#
# functions for displaying progress in parallel computing
# can be replaced with default progressmeter when https://github.com/timholy/ProgressMeter.jl/pull/157 is merged

import ProgressMeter: next!, finish!, cancel, update!
import Distributed: RemoteChannel, put!

struct ParallelProgress{C}
    channel::C
    n::Int
end

const PP_NEXT = -1
const PP_FINISH = -2
const PP_CANCEL = -3

next!(pp::ParallelProgress) = put!(pp.channel, PP_NEXT)
finish!(pp::ParallelProgress) = put!(pp.channel, PP_FINISH)
cancel(pp::ParallelProgress, args...; kw...) = put!(pp.channel, PP_CANCEL)
update!(pp::ParallelProgress, counter, color = nothing) = put!(pp.channel, counter)

function ParallelProgress(n::Int; kw...)
    channel = RemoteChannel(() -> Channel{Int}(n))
    progress = Progress(n; kw...)

    @async while progress.counter < progress.n
        f = take!(channel)
        if f == PP_NEXT
            next!(progress)
        elseif f == PP_FINISH
            finish!(progress)
            break
        elseif f == PP_CANCEL
            cancel(progress)
            break
        elseif f >= 0
            update!(progress, f)
        end
    end
    return ParallelProgress(channel, n)
end

struct MultipleChannel{C}
    channel::C
    id
end
put!(mc::MultipleChannel, x) = put!(mc.channel, (mc.id, x))

struct MultipleProgress{C}
    channel::C
    amount::Int
    lengths::Vector{Int}
end

Base.getindex(mp::MultipleProgress, n::Integer) = ParallelProgress(MultipleChannel(mp.channel, n), mp.lengths[n])
finish!(mp::MultipleProgress) = put!.([mp.channel], [(p, PP_FINISH) for p in 1:mp.amount])

function MultipleProgress(amount::Integer, lengths::Integer; kw...)
    MultipleProgress(amount, fill(lengths, amount); kw...)
end

function MultipleProgress(amount::Integer,
                          lengths::AbstractVector{<:Integer};
                          update_period = 0.1,
                          kws = [() for _ in 1:amount],
                          kw...)
    @assert amount == length(lengths) "`length(lengths)` must be equal to `amount`"

    total_length = sum(lengths)
    main_progress = Progress(total_length; offset=0, kw...)
    progresses = Union{Progress,Nothing}[nothing for _ in 1:amount]
    taken_offsets = Set(Int[])
    channel = RemoteChannel(() -> Channel{Tuple{Int,Int}}(max(2amount, 64)))

    max_offsets = 1

    # we must make sure that 2 progresses aren't updated at the same time
    @async begin
        while true

            (p, value) = take!(channel)

            # first time calling progress p
            if isnothing(progresses[p])
                # find first available offset
                offset = 1
                while offset in taken_offsets
                    offset += 1
                end
                max_offsets = max(max_offsets, offset)
                progresses[p] = Progress(lengths[p]; offset=offset, kw..., kws[p]...)
                push!(taken_offsets, offset)
            end

            if value == PP_NEXT
                next!(progresses[p])
                next!(main_progress)
            else
                prev_p_value = progresses[p].counter

                if value == PP_FINISH
                    finish!(progresses[p])
                elseif value == PP_CANCEL
                    cancel(progresses[p])
                elseif value >= 0
                    update!(progresses[p], value)
                end

                update!(main_progress,
                        main_progress.counter - prev_p_value + progresses[p].counter)
            end

            if progresses[p].counter >= lengths[p]
                delete!(taken_offsets, progresses[p].offset)
            end

            main_progress.counter >= total_length && break
        end

        print("\n" ^ max_offsets)
    end

    return MultipleProgress(channel, amount, lengths)
end


text:
using Distributed
addprocs(2)

@everywhere begin
    using MLJParticleSwarmOptimization
    using Random
    using Test
    using ComputationalResources
    using Distributions
    using EvoTrees
    using MLJBase
    using MLJTuning
    using StableRNGs
end

const PSO = MLJParticleSwarmOptimization

include("parameters.jl")
include("update.jl")
include("problems.jl")
include("strategies/basic.jl")
include("strategies/adaptive.jl")


text:
# Various Unicode functionality from the utf8proc library
module UTF8proc

import Base: show, showcompact, ==, string, symbol, isless

# also exported by Base:
export normalize_string, is_valid_char, is_assigned_char

# whether codepoints are valid Unicode
is_valid_char(c) = bool(ccall(:utf8proc_codepoint_valid, Cchar, (Int32,), c))

const UTF8PROC_NULLTERM  = (1<<0)
const UTF8PROC_STABLE    = (1<<1)
const UTF8PROC_COMPAT    = (1<<2)
const UTF8PROC_COMPOSE   = (1<<3)
const UTF8PROC_DECOMPOSE = (1<<4)
const UTF8PROC_IGNORE    = (1<<5)
const UTF8PROC_REJECTNA  = (1<<6)
const UTF8PROC_NLF2LS    = (1<<7)
const UTF8PROC_NLF2PS    = (1<<8)
const UTF8PROC_NLF2LF    = (UTF8PROC_NLF2LS | UTF8PROC_NLF2PS)
const UTF8PROC_STRIPCC   = (1<<9)
const UTF8PROC_CASEFOLD  = (1<<10)
const UTF8PROC_CHARBOUND = (1<<11)
const UTF8PROC_LUMP      = (1<<12)
const UTF8PROC_STRIPMARK = (1<<13)

let
    const p = Array(Ptr{Uint8}, 1)
    global utf8proc_map
    function utf8proc_map(s::String, flags::Integer)
        result = ccall(:utf8proc_map, Cssize_t,
                       (Ptr{Uint8}, Cssize_t, Ptr{Ptr{Uint8}}, Cint),
                       bytestring(s), 0, p, flags | UTF8PROC_NULLTERM)
        result < 0 && error(bytestring(ccall(:utf8proc_errmsg, Ptr{Uint8},
                                             (Cssize_t,), result)))
        a = ccall(:jl_ptr_to_array_1d, Vector{Uint8},
                  (Any, Ptr{Uint8}, Csize_t, Cint),
                  Vector{Uint8}, p[1], result, true)
        ccall(:jl_array_to_string, Any, (Any,), a)::ByteString
    end
end

function normalize_string(s::String; stable::Bool=false, compat::Bool=false, compose::Bool=true, decompose::Bool=false, stripignore::Bool=false, rejectna::Bool=false, newline2ls::Bool=false, newline2ps::Bool=false, newline2lf::Bool=false, stripcc::Bool=false, casefold::Bool=false, lump::Bool=false, stripmark::Bool=false)
    flags = 0
    stable && (flags = flags | UTF8PROC_STABLE)
    compat && (flags = flags | UTF8PROC_COMPAT)
    if decompose
        flags = flags | UTF8PROC_DECOMPOSE
    elseif compose
        flags = flags | UTF8PROC_COMPOSE
    elseif compat || stripmark
        throw(ArgumentError("compat=true or stripmark=true require compose=true or decompose=true"))
    end
    stripignore && (flags = flags | UTF8PROC_IGNORE)
    rejectna && (flags = flags | UTF8PROC_REJECTNA)
    newline2ls + newline2ps + newline2lf > 1 && throw(ArgumentError("only one newline conversion may be specified"))
    newline2ls && (flags = flags | UTF8PROC_NLF2LS)
    newline2ps && (flags = flags | UTF8PROC_NLF2PS)
    newline2lf && (flags = flags | UTF8PROC_NLF2LF)
    stripcc && (flags = flags | UTF8PROC_STRIPCC)
    casefold && (flags = flags | UTF8PROC_CASEFOLD)
    lump && (flags = flags | UTF8PROC_LUMP)
    stripmark && (flags = flags | UTF8PROC_STRIPMARK)
    utf8proc_map(s, flags)
end

function normalize_string(s::String, nf::Symbol)
    utf8proc_map(s, nf == :NFC ? (UTF8PROC_STABLE | UTF8PROC_COMPOSE) :
                    nf == :NFD ? (UTF8PROC_STABLE | UTF8PROC_DECOMPOSE) :
                    nf == :NFKC ? (UTF8PROC_STABLE | UTF8PROC_COMPOSE
                                   | UTF8PROC_COMPAT) :
                    nf == :NFKD ? (UTF8PROC_STABLE | UTF8PROC_DECOMPOSE
                                   | UTF8PROC_COMPAT) :
                    throw(ArgumentError(":$nf is not one of :NFC, :NFD, :NFKC, :NFKD")))
end

# returns UTF8PROC_CATEGORY code in 0..30 giving Unicode category
function category_code(c)
    # note: utf8proc returns 0, not UTF8PROC_CATEGORY_CN, for unassigned c
    c > 0x10FFFF && return 0x0000 # see utf8proc_get_property docs
    unsafe_load(ccall(:utf8proc_get_property, Ptr{Uint16}, (Int32,), c))
end

is_assigned_char(c) = category_code(c) != 0

# TODO: use UTF8PROC_CHARBOUND to extract graphemes from a string, e.g. to iterate over graphemes?

end # module

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia