text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################
# test/model.jl
#############################################################################

module TestModels

using JuMP
using Test

# Simple LP model not supporting Interval
MOIU.@model(
    SimpleLPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    ()
)

struct Optimizer
    a::Int
    b::Int
end
function opt_build(a::Int; b::Int = 1)
    return Optimizer(a, b)
end

# Custom set Nonnegative with bridge NonnegativeBridge
include("nonnegative_bridge.jl")

function _test_result_attributes(; test_empty = false)
    err = JuMP.OptimizeNotCalled()
    model = Model(() -> MOIU.MockOptimizer(SimpleLPModel{Float64}()))
    @variable(model, x)
    c = @constraint(model, x ≤ 0)
    @objective(model, Max, x)
    if test_empty
        optimize!(model)
        empty!(model)
    end
    @test_throws err JuMP.objective_value(model)
    @test_throws err JuMP.dual_objective_value(model)
    @test_throws err JuMP.objective_bound(model)
    @test_throws err JuMP.value(x)
    @test_throws err JuMP.value(c)
    @test_throws err JuMP.dual(c)
end

function test_result_attributes()
    return _test_result_attributes()
end

function test_result_attributes_after_empty()
    return _test_result_attributes(test_empty = true)
end

function fill_small_test_model!(model)
    # The model does not need to make sense, just use many different features.
    @variable(model, a[1:5] >= 0, Int)
    @variable(model, b[6:10], Bin)
    @variable(model, c[1:3] == 0)
    @variable(model, 10 <= d[1:3] <= 20)
    @constraint(model, con1, sum(a) + sum(b) <= 5)
    @constraint(model, con2, sum(b) >= 3)
    @constraint(model, con3, sum(d[1:2]) >= 5)
    @constraint(model, con4, sum(d) <= (sum(c) + 10))
    @objective(model, Max, sum(a) - sum(b) + sum(d))
    return model
end

function test_nooptimizer()
    err = NoOptimizer()
    model = Model()
    @variable(model, x)
    @test_throws err optimizer_index(x)
    cref = @constraint(model, x == 1)
    @test_throws err JuMP.optimizer_index(cref)
    @test_throws err JuMP.optimize!(model)
    @test_throws err JuMP.value(x)
    @test_throws err JuMP.value(cref)
    @test_throws err JuMP.dual(cref)
end

function test_empty!_model()
    model = Model()
    backend_type = typeof(backend(model))
    model.optimize_hook === nothing
    hook(m) = nothing
    JuMP.set_optimize_hook(model, hook)
    @test model.optimize_hook === hook
    @test fill_small_test_model!(model) === model
    @test_throws ErrorException fill_small_test_model!(model)
    @test empty!(model) === model
    @test model.optimize_hook === hook # empty! does not touch the hook
    @test isa(backend(model), backend_type)
    @test fill_small_test_model!(model) === model
end

function test_hygiene_variable()
    model_x = Model()
    @variable(model_x, x)
    model_y = Model()
    @variable(model_y, y)
    err = JuMP.VariableNotOwned(y)
    @test_throws err @constraint(model_x, y in MOI.EqualTo(1.0))
    @test_throws err @constraint(model_x, [x, y] in MOI.Zeros(2))
    @test_throws err @objective(model_x, Min, y)
end

function test_hygiene_linear()
    model_x = Model()
    @variable(model_x, x)
    model_y = Model()
    @variable(model_y, y)
    err = JuMP.VariableNotOwned(y)

    @test_throws err @constraint(model_x, x + y == 1)
    @test_throws err @constraint(model_x, [x, x + y] in MOI.Zeros(2))
    @test_throws err @constraint(model_x, [x + y, x] in MOI.Zeros(2))
    @test_throws err @objective(model_x, Min, x + y)
end

function test_hygiene_quadratic()
    model_x = Model()
    @variable(model_x, x)
    model_y = Model()
    @variable(model_y, y)
    err = JuMP.VariableNotOwned(y)

    @test_throws err @constraint(model_x, x * y >= 1)
    @test_throws err @constraint(model_x, [x, x * y] in MOI.Zeros(2))
    @test_throws err @constraint(model_x, [x * y, x] in MOI.Zeros(2))
    @test_throws err @constraint(model_x, x * x + x + y <= 1)
    @test_throws err @constraint(model_x, [x, x * x + x + y] in MOI.Zeros(2))
    @test_throws err @constraint(model_x, [x * x + x + y, x] in MOI.Zeros(2))
    @test_throws err @objective(model_x, Min, x * y)
    @test_throws err @objective(model_x, Min, x * x + x + y)
end

function test_hygiene_attribute()
    model_x = Model()
    @variable(model_x, x)
    model_y = Model()
    @variable(model_y, y)
    err = JuMP.VariableNotOwned(y)

    cy = @constraint(model_y, y in MOI.EqualTo(1.0))
    cerr = JuMP.ConstraintNotOwned(cy)

    @test_throws err MOI.get(model_x, MOI.VariablePrimalStart(), y)
    @test_throws cerr MOI.get(model_x, MOI.ConstraintPrimalStart(), cy)
    @test_throws err MOI.set(model_x, MOI.VariablePrimalStart(), y, 1.0)
    @test_throws cerr MOI.set(model_x, MOI.ConstraintPrimalStart(), cy, 1.0)
end

function test_optimize_hook()
    m = Model()
    @test m.optimize_hook === nothing
    called = false
    function hook(m)
        return called = true
    end
    JuMP.set_optimize_hook(m, hook)
    @test !called
    optimize!(m)
    @test called

    m = Model()
    err = ErrorException("Unrecognized keyword arguments: unexpected_arg")
    @test_throws err optimize!(m, unexpected_arg = 1)
    JuMP.set_optimize_hook(m, (m; my_new_arg = nothing) -> my_new_arg)
    @test optimize!(m) === nothing
    @test optimize!(m, my_new_arg = 1) == 1
end

function test_universal_fallback()
    m = Model()
    MOI.set(m, MOI.Test.UnknownModelAttribute(), 1)
    @test MOI.get(m, MOI.Test.UnknownModelAttribute()) == 1
end

function test_bridges_automatic()
    # optimizer not supporting Interval
    model = Model(() -> MOIU.MockOptimizer(SimpleLPModel{Float64}()))
    @test JuMP.bridge_constraints(model)
    @test JuMP.backend(model) isa MOIU.CachingOptimizer
    @test JuMP.backend(model).optimizer isa MOI.Bridges.LazyBridgeOptimizer
    @test JuMP.backend(model).optimizer.model isa MOIU.MockOptimizer
    @variable model x
    cref = @constraint model 0 <= x + 1 <= 1
    @test cref isa JuMP.ConstraintRef{
        JuMP.Model,
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.Interval{Float64},
        },
    }
    return JuMP.optimize!(model)
end

function test_bridges_automatic_with_cache()
    # Automatic bridging with cache for bridged model
    # optimizer not supporting Interval and not supporting `default_copy_to`
    model = Model(
        () -> MOIU.MockOptimizer(
            SimpleLPModel{Float64}(),
            needs_allocate_load = true,
        ),
    )
    @test JuMP.bridge_constraints(model)
    @test JuMP.backend(model) isa MOIU.CachingOptimizer
    @test JuMP.backend(model).optimizer isa MOI.Bridges.LazyBridgeOptimizer
    @test JuMP.backend(model).optimizer.model isa MOIU.CachingOptimizer
    @test JuMP.backend(model).optimizer.model.optimizer isa MOIU.MockOptimizer
    @variable model x
    err = ErrorException(
        "There is no `optimizer_index` as the optimizer is not " *
        "synchronized with the cached model. Call " *
        "`MOIU.attach_optimizer(model)` to synchronize it.",
    )
    @test_throws err optimizer_index(x)
    cref = @constraint model 0 <= x + 1 <= 1
    @test cref isa JuMP.ConstraintRef{
        JuMP.Model,
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.Interval{Float64},
        },
    }
    @test_throws err optimizer_index(cref)
    JuMP.optimize!(model)
    err = ErrorException(
        "There is no `optimizer_index` for $(typeof(index(cref))) " *
        "constraints because they are bridged.",
    )
    @test_throws err optimizer_index(cref)
end

function test_bridges_automatic_disabled()
    # Automatic bridging disabled with `bridge_constraints` keyword
    model = Model(
        () -> MOIU.MockOptimizer(SimpleLPModel{Float64}()),
        bridge_constraints = false,
    )
    @test !JuMP.bridge_constraints(model)
    @test JuMP.backend(model) isa MOIU.CachingOptimizer
    @test !(JuMP.backend(model).optimizer isa MOI.Bridges.LazyBridgeOptimizer)
    @variable model x
    err = ErrorException(
        "Constraints of type MathOptInterface.ScalarAffineFunction{Float64}-in-MathOptInterface.Interval{Float64} are not supported by the solver, try using `bridge_constraints=true` in the `JuMP.Model` constructor if you believe the constraint can be reformulated to constraints supported by the solver.",
    )
    @test_throws err @constraint model 0 <= x + 1 <= 1
end

function test_bridges_direct()
    # No bridge automatically added in Direct mode
    optimizer = MOIU.MockOptimizer(SimpleLPModel{Float64}())
    model = JuMP.direct_model(optimizer)
    @test !JuMP.bridge_constraints(model)
    @variable model x
    err = ErrorException(
        "Constraints of type MathOptInterface.ScalarAffineFunction{Float64}-in-MathOptInterface.Interval{Float64} are not supported by the solver.",
    )
    @test_throws err @constraint model 0 <= x + 1 <= 1
end

function mock_factory()
    mock = MOIU.MockOptimizer(
        MOIU.Model{Float64}(),
        eval_variable_constraint_dual = false,
    )
    function optimize!(mock)
        return MOIU.mock_optimize!(
            mock,
            [1.0],
            (MOI.SingleVariable, MOI.GreaterThan{Float64}) => [2.0],
        )
    end
    MOIU.set_mock_optimize!(mock, optimize!)
    return mock
end

function test_bridges_add_before_con_model_optimizer()
    model = Model(mock_factory)
    @variable(model, x)
    JuMP.add_bridge(model, NonnegativeBridge)
    c = @constraint(model, x in Nonnegative())
    JuMP.optimize!(model)
    @test 1.0 == @inferred JuMP.value(x)
    @test 1.0 == @inferred JuMP.value(c)
    @test 2.0 == @inferred JuMP.dual(c)
end

function test_bridges_add_before_con_set_optimizer()
    model = Model()
    @variable(model, x)
    c = @constraint(model, x in Nonnegative())
    JuMP.add_bridge(model, NonnegativeBridge)
    set_optimizer(model, mock_factory)
    JuMP.optimize!(model)
    @test 1.0 == @inferred JuMP.value(x)
    @test 1.0 == @inferred JuMP.value(c)
    @test 2.0 == @inferred JuMP.dual(c)
end

function test_bridges_add_after_con_model_optimizer()
    model = Model(mock_factory)
    @variable(model, x)
    flag = true
    try
        @constraint(model, x in Nonnegative())
        flag = false
    catch err
        @test err isa ErrorException
        # Rather than test a particular bridging error, just check that the
        # bridge explanation has been called. The sequence of errors could vary
        # between MOI versions.
        @test occursin("Nonnegative", err.msg)
        @test occursin("are not supported and cannot be bridged", err.msg)
    end
    @test flag
    JuMP.add_bridge(model, NonnegativeBridge)
    c = @constraint(model, x in Nonnegative())
    JuMP.optimize!(model)
    @test 1.0 == @inferred JuMP.value(x)
    @test 1.0 == @inferred JuMP.value(c)
    @test 2.0 == @inferred JuMP.dual(c)
end

function test_bridges_add_after_con_set_optimizer()
    err = MOI.UnsupportedConstraint{MOI.SingleVariable,Nonnegative}()
    model = Model()
    @variable(model, x)
    c = @constraint(model, x in Nonnegative())
    set_optimizer(model, mock_factory)
    @test_throws err JuMP.optimize!(model)
    JuMP.add_bridge(model, NonnegativeBridge)
    JuMP.optimize!(model)
    @test 1.0 == @inferred JuMP.value(x)
    @test 1.0 == @inferred JuMP.value(c)
    @test 2.0 == @inferred JuMP.dual(c)
end

function test_bridges_add_bridgeable_con_model_optimizer()
    model = Model(mock_factory)
    @variable(model, x)
    constraint = ScalarConstraint(x, Nonnegative())
    bc = BridgeableConstraint(constraint, NonnegativeBridge)
    c = add_constraint(model, bc)
    JuMP.optimize!(model)
    @test 1.0 == @inferred JuMP.value(x)
    @test 1.0 == @inferred JuMP.value(c)
    @test 2.0 == @inferred JuMP.dual(c)
end

function test_bridges_add_bridgeable_con_set_optimizer()
    model = Model()
    @variable(model, x)
    constraint = ScalarConstraint(x, Nonnegative())
    bc = BridgeableConstraint(constraint, NonnegativeBridge)
    c = add_constraint(model, bc)
    set_optimizer(model, mock_factory)
    JuMP.optimize!(model)
    @test 1.0 == @inferred JuMP.value(x)
    @test 1.0 == @inferred JuMP.value(c)
    @test 2.0 == @inferred JuMP.dual(c)
end

function test_bridge_graph_false()
    model = Model(mock_factory, bridge_constraints = false)
    @variable(model, x)
    @test_throws(
        ErrorException(
            "Cannot add bridge if `bridge_constraints` was set to `false` in " *
            "the `Model` constructor.",
        ),
        add_bridge(model, NonnegativeBridge)
    )
    @test_throws(
        ErrorException(
            "Cannot print bridge graph if `bridge_constraints` was set to " *
            "`false` in the `Model` constructor.",
        ),
        print_bridge_graph(model)
    )
    optimize!(model)
    @test 1.0 == @inferred value(x)
end

function test_bridge_graph_true()
    model = Model(mock_factory)
    @variable(model, x)
    add_bridge(model, NonnegativeBridge)
    @test sprint(print_bridge_graph, model) ==
          "Bridge graph with 0 variable nodes, 0 constraint nodes and 0 objective nodes.\n"
    c = @constraint(model, x in Nonnegative())
    @test sprint(print_bridge_graph, model) == replace(
        "Bridge graph with 1 variable nodes, 2 constraint nodes and 0 objective nodes.\n" *
        " [1] constrained variables in `$(Nonnegative)` are supported (distance 2) by adding free variables and then constrain them, see (1).\n" *
        " (1) `MOI.SingleVariable`-in-`$(Nonnegative)` constraints are bridged (distance 1) by $(NonnegativeBridge{Float64,MOI.SingleVariable}).\n" *
        " (2) `MOI.ScalarAffineFunction{Float64}`-in-`$(Nonnegative)` constraints are bridged (distance 1) by $(NonnegativeBridge{Float64,MOI.ScalarAffineFunction{Float64}}).\n",
        "MathOptInterface." => "MOI.",
    )
    optimize!(model)
    @test 1.0 == @inferred value(x)
    @test 1.0 == @inferred value(c)
    @test 2.0 == @inferred dual(c)
end

function test_solve_time()
    err = NoOptimizer()
    model = Model()
    @test_throws err solve_time(model)

    err = OptimizeNotCalled()
    model = Model(() -> MOIU.MockOptimizer(SimpleLPModel{Float64}()))
    @test_throws err solve_time(model)

    # TODO: Solved model
end

function test_solver_name()
    model = Model()
    @test "No optimizer attached." == @inferred JuMP.solver_name(model)

    model = Model(() -> MOIU.MockOptimizer(SimpleLPModel{Float64}()))
    @test "Mock" == @inferred JuMP.solver_name(model)
end

function test_set_silent()
    mock = MOIU.UniversalFallback(MOIU.Model{Float64}())
    model = Model(() -> MOIU.MockOptimizer(mock))
    @test JuMP.set_silent(model)
    @test MOI.get(backend(model), MOI.Silent())
    @test MOI.get(model, MOI.Silent())
    @test !JuMP.unset_silent(model)
    @test !MOI.get(backend(model), MOI.Silent())
    @test !MOI.get(model, MOI.Silent())
end

function test_set_optimizer_attribute()
    mock = MOIU.UniversalFallback(MOIU.Model{Float64}())
    model = Model(() -> MOIU.MockOptimizer(mock))
    @test JuMP.set_optimizer_attribute(model, "aaa", "bbb") == "bbb"
    @test MOI.get(backend(model), MOI.RawParameter("aaa")) == "bbb"
    @test MOI.get(model, MOI.RawParameter("aaa")) == "bbb"
end

function test_set_optimizer_attributes()
    mock = MOIU.UniversalFallback(MOIU.Model{Float64}())
    model = Model(() -> MOIU.MockOptimizer(mock))
    JuMP.set_optimizer_attributes(model, "aaa" => "bbb", "abc" => 10)
    @test MOI.get(model, MOI.RawParameter("aaa")) == "bbb"
    @test MOI.get(model, MOI.RawParameter("abc")) == 10
end

function test_get_optimizer_attribute()
    mock = MOIU.UniversalFallback(MOIU.Model{Float64}())
    model = Model(() -> MOIU.MockOptimizer(mock))
    @test JuMP.set_optimizer_attribute(model, "aaa", "bbb") == "bbb"
    @test JuMP.get_optimizer_attribute(model, "aaa") == "bbb"
end

function test_set_retrieve_time_limit()
    mock = MOIU.UniversalFallback(MOIU.Model{Float64}())
    model = Model(() -> MOIU.MockOptimizer(mock))
    JuMP.set_time_limit_sec(model, 12.0)
    @test JuMP.time_limit_sec(model) == 12.0
    JuMP.set_time_limit_sec(model, nothing)
    @test JuMP.time_limit_sec(model) === nothing
    JuMP.set_time_limit_sec(model, 12.0)
    JuMP.unset_time_limit_sec(model)
    @test JuMP.time_limit_sec(model) === nothing
end

struct DummyExtensionData
    model::JuMP.Model
end
function JuMP.copy_extension_data(
    data::DummyExtensionData,
    new_model::JuMP.AbstractModel,
    model::JuMP.AbstractModel,
)
    @test data.model === model
    return DummyExtensionData(new_model)
end
function dummy_optimizer_hook(::JuMP.AbstractModel) end

function copy_model_style_mode(use_copy_model, caching_mode, filter_mode)
    model = Model(caching_mode = caching_mode)
    model.optimize_hook = dummy_optimizer_hook
    data = DummyExtensionData(model)
    model.ext[:dummy] = data
    @variable(model, w[i = 1:2] ≥ 0)
    @variable(model, x ≥ 0, Bin)
    @variable(model, y ≤ 1, Int)
    @variable(model, z == 0)
    @constraint(model, cref, x + y == 1)
    @constraint(model, cref2[i = 1:2], w[i] + z == 1)

    if use_copy_model
        if filter_mode
            filter_constraints = (cr) -> cr != cref
            new_model, reference_map =
                JuMP.copy_model(model, filter_constraints = filter_constraints)
        else
            new_model, reference_map = JuMP.copy_model(model)
        end
    else
        new_model = copy(model)
        reference_map = Dict{
            Union{JuMP.VariableRef,JuMP.ConstraintRef},
            Union{JuMP.VariableRef,JuMP.ConstraintRef},
        }()
        reference_map[w[1]] = new_model[:w][1]
        reference_map[w[2]] = new_model[:w][2]
        reference_map[x] = new_model[:x]
        reference_map[y] = new_model[:y]
        reference_map[z] = new_model[:z]
        reference_map[cref] = new_model[:cref]
        reference_map[cref2[1]] = new_model[:cref2][1]
        reference_map[cref2[2]] = new_model[:cref2][2]
    end
    @test caching_mode == @inferred MOIU.mode(JuMP.backend(new_model))
    @test new_model.optimize_hook === dummy_optimizer_hook
    @test new_model.ext[:dummy].model === new_model
    x_new = reference_map[x]
    @test JuMP.owner_model(x_new) === new_model
    @test "x" == @inferred JuMP.name(x_new)
    y_new = reference_map[y]
    @test JuMP.owner_model(y_new) === new_model
    @test "y" == @inferred JuMP.name(y_new)
    z_new = reference_map[z]
    @test JuMP.owner_model(z_new) === new_model
    @test "z" == @inferred JuMP.name(z_new)
    if use_copy_model
        @test reference_map[JuMP.LowerBoundRef(x)] ==
              @inferred JuMP.LowerBoundRef(x_new)
        @test reference_map[JuMP.BinaryRef(x)] ==
              @inferred JuMP.BinaryRef(x_new)
        @test reference_map[JuMP.UpperBoundRef(y)] ==
              @inferred JuMP.UpperBoundRef(y_new)
        @test reference_map[JuMP.IntegerRef(y)] ==
              @inferred JuMP.IntegerRef(y_new)
        @test reference_map[JuMP.FixRef(z)] == @inferred JuMP.FixRef(z_new)
    end

    cref2_1_new = reference_map[cref2[1]]
    @test cref2_1_new.model === new_model
    @test "cref2[1]" == @inferred JuMP.name(cref2_1_new)
    cref2_2_new = reference_map[cref2[2]]
    @test cref2_2_new.model === new_model
    @test "cref2[2]" == @inferred JuMP.name(cref2_2_new)

    if filter_mode
        @test_throws KeyError JuMP.object_dictionary(new_model)[JuMP.name(cref)]
    else
        cref_new = reference_map[cref]
        @test cref_new.model === new_model
        @test "cref" == @inferred JuMP.name(cref_new)
    end
end

function test_copy_model_jump_auto()
    return copy_model_style_mode(true, MOIU.AUTOMATIC, false)
end

function test_compute_conflict()
    err = NoOptimizer()
    model = Model()
    @test_throws err compute_conflict!(model)
end

function test_copy_model_base_auto()
    return copy_model_style_mode(false, MOIU.AUTOMATIC, false)
end
function test_copy_model_jump_manual()
    return copy_model_style_mode(true, MOIU.MANUAL, false)
end
function test_copy_model_base_manual()
    return copy_model_style_mode(false, MOIU.MANUAL, false)
end

function test_copy_direct_mode()
    mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
    model = JuMP.direct_model(mock)
    @test_throws ErrorException JuMP.copy(model)
end

function test_copy_expr_aff()
    model = Model()
    @variable(model, x)
    @expression(model, ex, 2 * x + 1)
    new_model, ref_map = copy_model(model)
    @test ref_map[ex] == new_model[:ex]
    @test new_model[:ex] == 2 * new_model[:x] + 1
end

function test_copy_expr_quad()
    model = Model()
    @variable(model, x)
    @expression(model, ex, 2 * x^2 + x + 1)
    new_model, ref_map = copy_model(model)
    @test ref_map[ex] == new_model[:ex]
    @test new_model[:ex] == 2 * new_model[:x]^2 + new_model[:x] + 1
end

function test_copy_dense()
    model = Model()
    @variable(model, x[[:a, :b]])
    new_model, ref_map = copy_model(model)
    @test ref_map[x] == new_model[:x]
end

function test_copy_sparse()
    model = Model()
    @variable(model, x[i = 1:4; isodd(i)])
    new_model, ref_map = copy_model(model)
    @test ref_map[x] == new_model[:x]
end

function test_copy_object_fail()
    model = Model()
    @variable(model, x)
    model[:d] = Dict(x => 2)
    new_model, ref_map = copy_model(model)
    @test !haskey(new_model, :d)
end

function test_haskey()
    model = Model()
    @variable(model, p[i = 1:10] >= 0)
    @test haskey(model, :p)
    @test !haskey(model, :i)
end

function test_copy_refmap_expr()
    model = Model()
    @variable(model, x)
    @expression(model, expr[i = 1:2, j = 1:2], [i, j] * x)
    new_model, ref_map = copy_model(model)
    @test ref_map[expr] == new_model[:expr]
    for i in 1:2, j in 1:2
        @test new_model[:expr][i, j] == [i, j] * new_model[:x]
    end
end

function test_copy_dict_expr()
    model = Model()
    @variable(model, x)
    @expression(model, dictExpr[i = 1:2, j = 1:2], Dict(i => x, j => 2x))
    new_model, ref_map = copy_model(model)
    @test !haskey(new_model, :dictExpr)
end

function test_copy_filter()
    return copy_model_style_mode(true, MOIU.AUTOMATIC, true)
end

function test_copy_filter_array()
    model = Model()
    @variable(model, x[i = 1:2], container = Array)
    @constraint(model, cref[i = 1:2], x[i] == 1, container = Array)
    @test num_constraints(
        model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 2

    filter_constraints = (cr) -> cr != cref[1]
    new_model, reference_map =
        JuMP.copy_model(model, filter_constraints = filter_constraints)
    @test num_constraints(
        new_model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 1

    x1_new = reference_map[x[1]]
    @test JuMP.owner_model(x1_new) === new_model
    @test "x[1]" == @inferred JuMP.name(x1_new)

    cref_2_new = reference_map[cref[2]]
    @test cref_2_new.model === new_model
    @test "cref[2]" == @inferred JuMP.name(cref_2_new)
end

function test_copy_filter_denseaxisarray()
    model = Model()
    @variable(model, x[i = 1:2], container = DenseAxisArray)
    @constraint(model, cref[i = 1:2], x[i] == 1, container = DenseAxisArray)
    @test num_constraints(
        model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 2

    filter_constraints = (cr) -> cr != cref[1]
    new_model, reference_map =
        JuMP.copy_model(model, filter_constraints = filter_constraints)
    @test num_constraints(
        new_model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 1

    x1_new = reference_map[x[1]]
    @test JuMP.owner_model(x1_new) === new_model
    @test "x[1]" == @inferred JuMP.name(x1_new)

    cref_2_new = reference_map[cref[2]]
    @test cref_2_new.model === new_model
    @test "cref[2]" == @inferred JuMP.name(cref_2_new)
end

function test_copy_filter_sparseaxisarray()
    model = Model()
    @variable(model, x[i = 1:2], container = SparseAxisArray)
    @constraint(model, cref[i = 1:2], x[i] == 1, container = SparseAxisArray)
    @test num_constraints(
        model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 2

    filter_constraints = (cr) -> cr != cref[1]
    new_model, reference_map =
        JuMP.copy_model(model, filter_constraints = filter_constraints)
    @test num_constraints(
        new_model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 1

    x1_new = reference_map[x[1]]
    @test JuMP.owner_model(x1_new) === new_model
    @test "x[1]" == @inferred JuMP.name(x1_new)

    cref_2_new = reference_map[cref[2]]
    @test cref_2_new.model === new_model
    @test "cref[2]" == @inferred JuMP.name(cref_2_new)
end

function test_copy_conflict()
    model = Model()
    @variable(model, x[i = 1:2], container = SparseAxisArray)
    @constraint(model, cref[i = 1:2], x[i] == 1, container = SparseAxisArray)
    @test num_constraints(
        model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 2

    set_optimizer(
        model,
        () -> MOIU.MockOptimizer(
            MOIU.Model{Float64}(),
            eval_objective_value = false,
        ),
    )
    JuMP.optimize!(model)

    mockoptimizer = JuMP.backend(model).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.INFEASIBLE)
    MOI.set(mockoptimizer, MOI.ConflictStatus(), MOI.CONFLICT_FOUND)
    MOI.set(
        mockoptimizer,
        MOI.ConstraintConflictStatus(),
        JuMP.optimizer_index(cref[1]),
        MOI.IN_CONFLICT,
    )
    MOI.set(
        mockoptimizer,
        MOI.ConstraintConflictStatus(),
        JuMP.optimizer_index(cref[2]),
        MOI.NOT_IN_CONFLICT,
    )

    new_model, reference_map = JuMP.copy_conflict(model)
    @test num_constraints(
        new_model,
        GenericAffExpr{Float64,VariableRef},
        MOI.EqualTo{Float64},
    ) == 1

    x1_new = reference_map[x[1]]
    @test JuMP.owner_model(x1_new) === new_model
    @test "x[1]" == @inferred JuMP.name(x1_new)

    cref_1_new = reference_map[cref[1]]
    @test cref_1_new.model === new_model
    @test "cref[1]" == @inferred JuMP.name(cref_1_new)
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        @testset "$(name)" begin
            getfield(@__MODULE__, name)()
        end
    end
end

end  # module TestModels

TestModels.runtests()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: dlm
using AWS.Compat
using AWS.UUIDs

"""
    create_lifecycle_policy(description, execution_role_arn, policy_details, state)
    create_lifecycle_policy(description, execution_role_arn, policy_details, state, params::Dict{String,<:Any})

Creates a policy to manage the lifecycle of the specified AWS resources. You can create up
to 100 lifecycle policies.

# Arguments
- `description`: A description of the lifecycle policy. The characters ^[0-9A-Za-z _-]+ are
  supported.
- `execution_role_arn`: The Amazon Resource Name (ARN) of the IAM role used to run the
  operations specified by the lifecycle policy.
- `policy_details`: The configuration details of the lifecycle policy.
- `state`: The desired activation state of the lifecycle policy after creation.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Tags"`: The tags to apply to the lifecycle policy during creation.
"""
create_lifecycle_policy(Description, ExecutionRoleArn, PolicyDetails, State; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("POST", "/policies", Dict{String, Any}("Description"=>Description, "ExecutionRoleArn"=>ExecutionRoleArn, "PolicyDetails"=>PolicyDetails, "State"=>State); aws_config=aws_config)
create_lifecycle_policy(Description, ExecutionRoleArn, PolicyDetails, State, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("POST", "/policies", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Description"=>Description, "ExecutionRoleArn"=>ExecutionRoleArn, "PolicyDetails"=>PolicyDetails, "State"=>State), params)); aws_config=aws_config)

"""
    delete_lifecycle_policy(policy_id)
    delete_lifecycle_policy(policy_id, params::Dict{String,<:Any})

Deletes the specified lifecycle policy and halts the automated operations that the policy
specified.

# Arguments
- `policy_id`: The identifier of the lifecycle policy.

"""
delete_lifecycle_policy(policyId; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("DELETE", "/policies/$(policyId)/"; aws_config=aws_config)
delete_lifecycle_policy(policyId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("DELETE", "/policies/$(policyId)/", params; aws_config=aws_config)

"""
    get_lifecycle_policies()
    get_lifecycle_policies(params::Dict{String,<:Any})

Gets summary information about all or the specified data lifecycle policies. To get
complete information about a policy, use GetLifecyclePolicy.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"policyIds"`: The identifiers of the data lifecycle policies.
- `"resourceTypes"`: The resource type.
- `"state"`: The activation state.
- `"tagsToAdd"`: The tags to add to objects created by the policy. Tags are strings in the
  format key=value. These user-defined tags are added in addition to the AWS-added lifecycle
  tags.
- `"targetTags"`: The target tag for a policy. Tags are strings in the format key=value.
"""
get_lifecycle_policies(; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("GET", "/policies"; aws_config=aws_config)
get_lifecycle_policies(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("GET", "/policies", params; aws_config=aws_config)

"""
    get_lifecycle_policy(policy_id)
    get_lifecycle_policy(policy_id, params::Dict{String,<:Any})

Gets detailed information about the specified lifecycle policy.

# Arguments
- `policy_id`: The identifier of the lifecycle policy.

"""
get_lifecycle_policy(policyId; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("GET", "/policies/$(policyId)/"; aws_config=aws_config)
get_lifecycle_policy(policyId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("GET", "/policies/$(policyId)/", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags for the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    tag_resource(tags, resource_arn)
    tag_resource(tags, resource_arn, params::Dict{String,<:Any})

Adds the specified tags to the specified resource.

# Arguments
- `tags`: One or more tags.
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

"""
tag_resource(Tags, resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("POST", "/tags/$(resourceArn)", Dict{String, Any}("Tags"=>Tags); aws_config=aws_config)
tag_resource(Tags, resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Tags"=>Tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes the specified tags from the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.
- `tag_keys`: The tag keys.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_lifecycle_policy(policy_id)
    update_lifecycle_policy(policy_id, params::Dict{String,<:Any})

Updates the specified lifecycle policy.

# Arguments
- `policy_id`: The identifier of the lifecycle policy.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: A description of the lifecycle policy.
- `"ExecutionRoleArn"`: The Amazon Resource Name (ARN) of the IAM role used to run the
  operations specified by the lifecycle policy.
- `"PolicyDetails"`: The configuration of the lifecycle policy. You cannot update the
  policy type or the resource type.
- `"State"`: The desired activation state of the lifecycle policy after creation.
"""
update_lifecycle_policy(policyId; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("PATCH", "/policies/$(policyId)"; aws_config=aws_config)
update_lifecycle_policy(policyId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = dlm("PATCH", "/policies/$(policyId)", params; aws_config=aws_config)

ext: jl


text:
module MapLight

import Base.show

using JSON
using HttpCommon
using Requests

export bill_positions,
       bill_search,
       bill_list,
       organization_positions,
       organization_search,
       authenticate

include("error.jl")
include("endpoint.jl")
include("auth.jl")
include("api.jl")

end


text:
module TablesDemo

using Reexport
using Compat
@reexport using AbstractTables
@reexport using DataStreams
@reexport using CSV
# @reexport using NullableArrays

export  Table

# Table
include("table/utils.jl")
include("table/typedef.jl")
include("table/traits.jl")
include("table/constructors.jl")
include("table/primitives.jl")
include("table/indexing.jl")
include("table/io.jl")

# query
include("query/interface.jl")

end # module TablesDemo


text:
using Test
using NSGAII

include("./cornercases.jl")

@testset "Domination" begin
    L=[[2,1,1],[1,2,2],[1,0,0]]
    fitness(x::FloatIndividual)=0
    cfg=Cambrian.get_config("./test.yaml";n_population= 3,n_offsprings=5,d_fitness=3)
    e=NSGA2Evolution{FloatIndividual}(cfg,fitness)
    for i in 1:3
        copyto!(e.population[i].fitness,L[i])
    end
    @test dominates(e,e.population[1],e.population[2])==false
    @test dominates(e,e.population[2],e.population[1])==false
    @test dominates(e,e.population[1],e.population[3])==true
    @test dominates(e,e.population[3],e.population[1])==false
end

@testset "fastNonDominatedSort" begin
    fitness1,ranks1=pop1D()
    fitness2,ranks2=pop3D()
    fitness(x::FloatIndividual)=0
    cfg=Cambrian.get_config("./test.yaml";n_population= 20,n_offsprings=5,d_fitness=1)
    e=NSGA2Evolution{FloatIndividual}(cfg,fitness)
    for i in 1:20
        copyto!(e.population[i].fitness,fitness1[i])
    end
    oldPop1D=copy(e.population)
    fastNonDominatedSort!(e)
    @test length(e.population)==length(oldPop1D)
    for x in oldPop1D
        @test x in e.population
    end
    for i in 1:e.config.n_population
        @test ranks1[i]==e.rank[objectid(e.population[i])]
    end

    cfg=Cambrian.get_config("./test.yaml";n_population= 16,n_offsprings=5,d_fitness=3)
    e=NSGA2Evolution{FloatIndividual}(cfg,fitness)
    for i in 1:16
        copyto!(e.population[i].fitness,fitness2[i])
    end
    oldPop3D=copy(e.population)
    fastNonDominatedSort!(e)
    @test length(e.population)==length(oldPop3D)
    for x in oldPop3D
        @test x in e.population
    end
    for i in 1:e.config.n_population
        @test ranks2[i]==e.rank[objectid(e.population[i])]
    end
end


text:
ind_others = [1,2]
h = 4/2
w = 2/2
r = w

ego_x,ego_y,ego_θ = 0, 0, 0
oth_x,oth_y,oth_θ = 6, 0, 0

# oth_xs = zeros(length(ind_others))
# oth_ys = zeros(length(ind_others))
# oth_θs = zeros(length(ind_others))
# # for (i,ind) in enumerate(ind_others)
# #     oth_xs[i],oth_ys[i],oth_θs[i] = scene[ind].state.posG.x, scene[ind].state.posG.y, scene[ind].state.posG.θ
# # end
#
# # TEMP
# oth_xs = [0 6 6 0];
# oth_ys = [3 0 3 -3];
# oth_θs = [0 0 0 0];

# @time begin
min_dis = Inf
for i in [-1,0,1]
    global min_dis
    for j in [-1,0,1]
        dis = sqrt(
                  ( (ego_x + i*r*cos(ego_θ)) - (oth_x+ j*r*cos(oth_θ)) )^2
                 +( (ego_y + i*r*sin(ego_θ)) - (oth_y+ j*r*sin(oth_θ)) )^2
                 )     - 2*r
        min_dis = min(dis, min_dis)
    end
end
# end

println("min distance: ", min_dis)


text:
# ProgressPrinter

"""
    ProgressPrinter()

Prints a progress bar of the currently running epoch.
"""
mutable struct ProgressPrinter <: Callback
    p::Union{Nothing,Progress}
end
ProgressPrinter() = ProgressPrinter(nothing)
Base.show(io::IO, ::ProgressPrinter) = print(io, "ProgressPrinter()")

function on(::EpochBegin,
        phase::Phase,
        cb::ProgressPrinter,
        learner)
    e = learner.cbstate.history.epochs + 1
    cb.p = Progress(numsteps(learner, phase), "Epoch $(e) $(phase): ")
end

on(::BatchEnd, ::Phase, cb::ProgressPrinter, learner) = next!(cb.p)

runafter(::ProgressPrinter) = (Recorder,)
stateaccess(::ProgressPrinter) = (data = Read(), cbstate = (history = Read()),)

"""
    MetricsPrinter()

Prints metrics after every epoch. Relies on [`Metrics`](#).
"""
struct MetricsPrinter <: Callback end

function on(::EpochEnd,
        phase::Phase,
        cb::MetricsPrinter,
        learner)
    mvhistory = learner.cbstate.metricsepoch[phase]
    epoch = learner.cbstate.history.epochs
    print_epoch_table(mvhistory, epoch, phase)
end

function print_epoch_table(mvhistory, epoch, phase)
    header = vcat(["Phase", "Epoch"], string.(keys(mvhistory)))
    vals = [last(mvhistory, key) |> last for key in keys(mvhistory)]
    data = reshape(vcat([string(typeof(phase)), epoch], vals), 1, :)
    pretty_table(data, header, formatters = PrettyTables.ft_round(5))
end

stateaccess(::MetricsPrinter) = (; cbstate = (metricsepoch = Read(), history = Read()))
runafter(::MetricsPrinter) = (Metrics,)

# StopOnNaNLoss

"""
    StopOnNaNLoss()

Stops the training when a NaN loss is encountered.
"""
struct StopOnNaNLoss <: Callback end

function on(::BackwardEnd, ::AbstractTrainingPhase, ::StopOnNaNLoss, learner)
    !isnan(learner.batch.loss) || throw(CancelFittingException("Encountered NaN loss"))
end

stateaccess(::StopOnNaNLoss) = (batch = (loss = Read()),)

# Early stopping
"""
    EarlyStopping(patience)

Stops training if validation loss hasn't decreased in `patience` epochs.
"""
mutable struct EarlyStopping <: Callback
    patience::Int
    waited::Int
    lowest::Float64
end
EarlyStopping(patience) = EarlyStopping(patience, 0, Inf64)

function on(::EpochEnd, ::ValidationPhase, cb::EarlyStopping, learner)
    valloss = last(learner.cbstate.metricsepoch[phase], :Loss)[2]
    if (valloss > cb.lowest)
        if !(cb.waited < cb.patience)
            throw(CancelFittingException("Validation loss did not improve for $(cb.patience) epochs"))
        else
            cb.waited += 1
        end
    else
        cb.waited = 0
        cb.lowest = valloss
    end
end

stateaccess(::EarlyStopping) = (callbacks = Read(), cbstate = (; metricsepoch = Read()))

"""
    ToGPU()

Callback that moves model and batch data to the GPU during training.
"""
struct ToGPU <: Callback end

function on(::EpochBegin, ::Phase, ::ToGPU, learner)
    model!(learner, gpu(learner.model))
end

stateaccess(::ToGPU) = (
    model = Write(),
    params = Write(),
    batch = (xs = Write(), ys = Write()),
)

function on(::BatchBegin, ::Phase, cb::ToGPU, learner)
    learner.batch.xs = gpu(learner.batch.xs)
    learner.batch.ys = gpu(learner.batch.ys)
end

function garbagecollect()
    GC.gc()
    if Base.Sys.islinux()
        ccall(:malloc_trim, Cvoid, (Cint,), 0)
    end
end

"""
    GarbageCollect(nsteps)

Every `nsteps` steps, forces garbage collection.
Use this if you get memory leaks from, for example,
parallel data loading.

Performs an additional C-call on Linux systems that can
sometimes help.
"""
function GarbageCollect(nsteps::Int = 100)
    return throttle(
        CustomCallback((learner) -> garbagecollect(), BatchEnd, Phase),
        BatchEnd(),
        freq = nsteps)
end


text:
# --------------------- #
# Post-estimation tools #
# --------------------- #

function optimal_W(mf::Function, theta::Vector, k::RobustVariance)
    h = mf(theta)
    n = size(h, 1)
    S = vcov(h, k) * n
    W = pinv(S)
    W
end

status(me::MomentEstimator) = me.r.status
StatsBase.coef(me::MomentEstimator) = me.r.coef
objval(me::MomentEstimator) = me.r.objval
momentfunction(me::MomentEstimator) = me.e.mf(coef(me))
jacobian(me::MomentEstimator) = me.e.Dmf(coef(me))
mfvcov(me::MomentEstimator, k::RobustVariance) = vcov(momentfunction(me), k)
nobs(me::MomentEstimator) = me.r.nobs
npar(me::MomentEstimator) = me.r.npar
nmom(me::MomentEstimator) = me.r.nmom
df(me::MomentEstimator) = nmom(me) - npar(me)
z_stats(me::MomentEstimator, k::RobustVariance) = coef(me) ./ stderr(me, k)
p_values(me::MomentEstimator, k::RobustVariance) = 2*ccdf(Normal(), z_stats(me, k))
shat(me::GMMEstimator, k::RobustVariance) = mfvcov(me, k)
optimal_W(me::GMMEstimator, k::RobustVariance) = pinv(full(shat(me, k)*nobs(me)))

StatsBase.vcov(me::MomentEstimator, k::RobustVariance) = vcov(me, k, me.e.mgr)
StatsBase.vcov(me::MomentEstimator) = vcov(me, me.e.mgr.k, me.e.mgr)

function StatsBase.vcov(me::MomentEstimator, k::RobustVariance, mgr::TwoStepGMM)
    G = jacobian(me)
    n = nobs(me)
    p = npar(me)
    S = shat(me, k)
    ## Use the general form of the variance covariance matrix
    ## that gives the correct covariance even when S \not Var(\sqrt{N}
    ## A = pinv(G'*pinv(S)*G)
    ## B = G'*pinv(S)**G
    (n.^2/(n-p))*pinv(G'*pinv(S)*G)
end

function StatsBase.vcov(me::MomentEstimator, k::RobustVariance, mgr::OneStepGMM)
    G = jacobian(me)
    n = nobs(me)
    p = npar(me)
    S = shat(me, k)
    W = me.e.W
    ## Use the general form of the variance covariance matrix
    ## that gives the correct covariance even when S \not Var(\sqrt{N}
    A = pinv(G'*W*G)
    B = G'*W*S*W*G
    (n.^2/(n-p))*A*B*A
end

function StatsBase.stderr(me::MomentEstimator)
    sqrt(diag(vcov(me, me.e.mgr.k, me.e.mgr)))
end

function StatsBase.stderr(me::MomentEstimator, mgr::IterationManager)
    sqrt(diag(vcov(me, mgr.k, mgr)))
end

function StatsBase.stderr(me::MomentEstimator, k::RobustVariance)
    sqrt(diag(vcov(me, k, me.e.mgr)))
end

function J_test(me::GMMEstimator, k::RobustVariance=me.e.mgr.k)
    g = mean(momentfunction(me), 1)
    S = pinv(shat(me, k))
    j = (nobs(me)*(g*S*g'))[1]
    p = df(me) > 0 ? ccdf(Chisq(df(me)), j) : NaN
    # sometimes p is garbage, so we clamp it to be within reason
    return j, clamp(p, eps(), Inf)
end

function StatsBase.coeftable(me::MomentEstimator,
                             k::RobustVariance=me.e.mgr.k)
    cc = coef(me)
    se = stderr(me, k)
    zz = z_stats(me, k)
    CoefTable(hcat(cc, se, zz, 2.0*ccdf(Normal(), abs(zz))),
              ["Estimate", "Std.Error", "z value", "Pr(>|z|)"],
              ["x$i" for i = 1:npar(me)],
              4)
end


text:
@testset "Conversion" begin
    P, Q = Chebyshev, ChebyshevU
    x = variable()

    for _ in 1:10
        p = P(rand(1:5, 5))
        q = P(rand(1:5, 5))
        @test convert(P, convert(Q, p))(x) ≈ p(x)
        @test convert(Q, convert(P, q))(x) ≈ q(x)
    end
end

@testset "Evaluation" begin
    p0 = ChebyshevU([1])  # 1
    p1 = ChebyshevU([0, 1]) # 2x
    p2 = ChebyshevU([0, 0, 1]) # 4x^2 - 1
    p3 = ChebyshevU([0, 0, 0, 1]) # 8x^3 - 4x
    p4 = ChebyshevU([0, 0, 0, 0, 1]) # 16x^4 - 12x^2 + 1

    x = 1 / 2
    @test p0(x) ≈ 1
    @test p1(x) ≈ 2x
    @test p2(x) ≈ 4x^2 - 1
    @test p3(x) ≈ 8x^3 - 4x
    @test p4(x) ≈ 16x^4 - 12x^2 + 1
    @test ChebyshevU([1, 2, 3, 4])(x) ≈ 1p0(x) + 2p1(x) + 3p2(x) + 4p3(x)
end

@testset "roots" begin
    for n in 4:2:8
        p = ChebyshevU(vcat(zeros(n), 1))
        @test cos.((n:-1:1) / (n + 1) * pi) ≈ roots(p)
    end
end


text:
@testset "Rectifier square" begin

r = Rectifier("squared")
@test r.T == "squared"

    x = 2.5
@test abs(r(x) - x^2)<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), x) - grad_x(r, x) ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - hess_x(r, x) ) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), x) - grad_x_logeval(r, x) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [x])[1,1] - hess_x_logeval(r, x) ) < 1e-10

x = -2.5
@test abs(r(x) - x^2)<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), x) - grad_x(r, x) ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - hess_x(r, x) ) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), x) - grad_x_logeval(r, x) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [x])[1,1] - hess_x_logeval(r, x) ) < 1e-10

end

@testset "Rectifier exponential" begin

r = Rectifier("exponential")
@test r.T == "exponential"

    x = 2.46
@test abs(r(x) - exp(x))<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), x) - grad_x(r, x) ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - hess_x(r, x) ) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), x) - grad_x_logeval(r, x) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [x])[1,1] - hess_x_logeval(r, x) ) < 1e-10

x = -2.46
@test abs(r(x) - exp(x))<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), x) - grad_x(r, x) ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - hess_x(r, x) ) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), x) - grad_x_logeval(r, x) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [x])[1,1] - hess_x_logeval(r, x) ) < 1e-10

end

@testset "Rectifier softplus" begin

r = Rectifier("softplus")
@test r.T == "softplus"

    x = 2.5
    a = log(2)
@test abs(r(x) - (log(1 + exp(-abs(a*x))) + max(a*x, 0))/a)<1e-10

#Test inversion
@test abs(inverse(r,r(x))-x)<1e-10
@test abs(r(inverse(r,x))-x)<1e-10

#Test inversion with x large
xlarge = 10.0^6
@test abs(inverse(r,xlarge) -xlarge)<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), x) - grad_x(r, x) ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - hess_x(r, x) ) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), x) - grad_x_logeval(r, x) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [x])[1,1] - hess_x_logeval(r, x) ) < 1e-10

x = -2.5
a = log(2)
@test abs(r(x) - (log(1 + exp(-abs(a*x))) + max(a*x, 0))/a)<1e-10

#Test inversion with x large
xlarge = 10.0^6
@test abs(inverse(r,xlarge) -xlarge)<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), x) - grad_x(r, x) ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - hess_x(r, x) ) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), x) - grad_x_logeval(r, x) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [x])[1,1] - hess_x_logeval(r, x) ) < 1e-10

end

@testset "Rectifier explinearunit" begin

r = Rectifier("explinearunit")
@test r.T == "explinearunit"

## Test x<0
xneg = -0.5
@test abs(r(xneg) - exp(xneg))<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), xneg) - grad_x(r, xneg) ) < 1e-10
@test abs(ForwardDiff.derivative(y->r(y), xneg) - exp(xneg) ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),xneg) - hess_x(r, xneg) ) < 1e-10
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),xneg) - exp(xneg) ) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), xneg) - grad_x_logeval(r, xneg) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [xneg])[1,1] - hess_x_logeval(r, xneg) ) < 1e-10

## Test x>=0
x = 2.5
@test abs(r(x) - (x + 1.0))<1e-10

#Test inversion
@test abs(inverse(r,r(x))-x)<1e-10
@test abs(r(inverse(r,x))-x)<1e-10

# Test gradient
@test abs(ForwardDiff.derivative(y->r(y), x) - grad_x(r, x) ) < 1e-10
@test abs(ForwardDiff.derivative(y->r(y), x) - 1.0 ) < 1e-10

# Test hessian
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - hess_x(r, x) ) < 1e-10
@test abs(ForwardDiff.derivative(z->ForwardDiff.derivative(y->r(y), z),x) - 0.0) < 1e-10

# Test gradient of log evaluation
@test abs(ForwardDiff.derivative(y->log(r(y)), x) - grad_x_logeval(r, x) ) < 1e-10

# Test hessian of log evaluation
@test abs(ForwardDiff.hessian(y->log(r(y[1])), [x])[1,1] - hess_x_logeval(r, x) ) < 1e-10

end


text:
using BinDeps

@BinDeps.setup

libstinger_core = library_dependency("libstinger_core")
libstinger_alg = library_dependency("libstinger_alg")

provides(Sources,
    URI("https://github.com/rohitvarkey/stinger/archive/stingergraphs.jl_v0.0.1.tar.gz"),
    [libstinger_core, libstinger_alg], unpacked_dir = "stinger-stingergraphs.jl_v0.0.1")

stingerbuilddir = joinpath(BinDeps.depsdir(libstinger_core), "src", "stinger-stingergraphs.jl_v0.0.1", "build")
prefix = joinpath(BinDeps.depsdir(libstinger_core), "usr")
provides(BuildProcess,
    (@build_steps begin
        GetSources(libstinger_core)
        CreateDirectory(stingerbuilddir)
        CreateDirectory(joinpath(prefix, "lib"))
        @build_steps begin
            ChangeDirectory(stingerbuilddir)
            `cmake .. -DBUILD_TESTING=OFF`
            `make`
            FileRule(joinpath(prefix, "lib" , "libstinger_core.$(Base.Libdl.dlext)"), @build_steps begin
                `cp lib/libstinger_core.$(Base.Libdl.dlext) $prefix/lib/`
            end)
            FileRule(joinpath(prefix, "lib" , "libstinger_alg.$(Base.Libdl.dlext)"), @build_steps begin
                `cp lib/libstinger_alg.$(Base.Libdl.dlext) $prefix/lib/`
            end)
        end
    end), [libstinger_core, libstinger_alg])

@BinDeps.install Dict(
    :libstinger_core => :libstinger_core,
    :libstinger_alg => :libstinger_alg)


text:
###############################################################################
#
#   fmpz_poly.jl : Flint polynomials over fmpz
#
###############################################################################

export FmpzPolyRing, fmpz_poly, cyclotomic, theta_qexp, eta_qexp, cos_minpoly,
       swinnerton_dyer, signature

###############################################################################
#
#   Data type and parent methods
#
###############################################################################

parent_type(::Type{fmpz_poly}) = FmpzPolyRing

elem_type(::Type{FmpzPolyRing}) = fmpz_poly

base_ring(a::FmpzPolyRing) = a.base_ring

parent(a::fmpz_poly) = a.parent

var(a::FmpzPolyRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

length(x::fmpz_poly) = ccall((:fmpz_poly_length, :libflint), Int,
                             (Ref{fmpz_poly},), x)

function coeff(x::fmpz_poly, n::Int)
   n < 0 && throw(DomainError())
   z = fmpz()
   ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
               (Ref{fmpz}, Ref{fmpz_poly}, Int), z, x, n)
   return z
end

zero(a::FmpzPolyRing) = a(0)

one(a::FmpzPolyRing) = a(1)

gen(a::FmpzPolyRing) = a([zero(base_ring(a)), one(base_ring(a))])

isgen(x::fmpz_poly) = ccall((:fmpz_poly_is_x, :libflint), Bool,
                            (Ref{fmpz_poly},), x)

function deepcopy_internal(a::fmpz_poly, dict::ObjectIdDict)
   z = fmpz_poly(a)
   z.parent = parent(a)
   return z
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(a::fmpz_poly) = canonical_unit(lead(a))

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::fmpz_poly)
   if length(x) == 0
      print(io, "0")
   else
      cstr = ccall((:fmpz_poly_get_str_pretty, :libflint), Ref{UInt8},
          (Ref{fmpz_poly}, Ref{UInt8}), x, string(var(parent(x))))

      print(io, unsafe_string(cstr))

      ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
   end
end

function show(io::IO, p::FmpzPolyRing)
   print(io, "Univariate Polynomial Ring in ")
   print(io, string(var(p)))
   print(io, " over ")
   show(io, p.base_ring)
end

show_minus_one(::Type{fmpz_poly}) = show_minus_one(fmpz)

###############################################################################
#
#   Unary operations
#
###############################################################################

function -(x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_neg, :libflint), Void,
         (Ref{fmpz_poly}, Ref{fmpz_poly}), z, x)
   return z
end

###############################################################################
#
#   Binary operations
#
###############################################################################

function +(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpz_poly_add, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly},  Ref{fmpz_poly}),
               z, x, y)
   return z
end

function -(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpz_poly_sub, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly},  Ref{fmpz_poly}),
               z, x, y)
   return z
end

function *(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpz_poly_mul, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly},  Ref{fmpz_poly}),
               z, x, y)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fmpz_poly)
   z = parent(y)()
   ccall((:fmpz_poly_scalar_mul_si, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, y, x)
   return z
end

function *(x::fmpz, y::fmpz_poly)
   z = parent(y)()
   ccall((:fmpz_poly_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz}), z, y, x)
   return z
end

function +(x::fmpz_poly, y::Int)
   z = parent(x)()
   ccall((:fmpz_poly_add_si, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, x, y)
   return z
end

function +(x::fmpz_poly, y::fmpz)
   z = parent(x)()
   ccall((:fmpz_poly_add_fmpz, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz}), z, x, y)
   return z
end

function -(x::fmpz_poly, y::Int)
   z = parent(x)()
   ccall((:fmpz_poly_sub_si, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, x, y)
   return z
end

function -(x::fmpz_poly, y::fmpz)
   z = parent(x)()
   ccall((:fmpz_poly_sub_fmpz, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz}), z, x, y)
   return z
end

function -(x::Int, y::fmpz_poly)
   z = parent(y)()
   ccall((:fmpz_poly_si_sub, :libflint), Void,
                (Ref{fmpz_poly}, Int, Ref{fmpz_poly}), z, x, y)
   return z
end

function -(x::fmpz, y::fmpz_poly)
   z = parent(y)()
   ccall((:fmpz_poly_fmpz_sub, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz}, Ref{fmpz_poly}), z, x, y)
   return z
end

+(x::Int, y::fmpz_poly) = y + x

+(x::fmpz, y::fmpz_poly) = y + x

*(x::fmpz_poly, y::Int) = y*x

*(x::fmpz_poly, y::fmpz) = y*x

+(x::Integer, y::fmpz_poly) = y + fmpz(x)

-(x::Integer, y::fmpz_poly) = fmpz(x) - y

*(x::Integer, y::fmpz_poly) = fmpz(x)*y

+(x::fmpz_poly, y::Integer) = x + fmpz(y)

-(x::fmpz_poly, y::Integer) = x - fmpz(y)

*(x::fmpz_poly, y::Integer) = fmpz(y)*x

###############################################################################
#
#   Powering
#
###############################################################################

function ^(x::fmpz_poly, y::Int)
   y < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpz_poly_pow, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Int),
               z, x, y)
   return z
end

###############################################################################
#
#   Comparisons
#
###############################################################################

function ==(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   return ccall((:fmpz_poly_equal, :libflint), Bool,
                (Ref{fmpz_poly}, Ref{fmpz_poly}), x, y)
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fmpz_poly, y::fmpz)
   if length(x) > 1
      return false
   elseif length(x) == 1
      z = fmpz()
      ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
                       (Ref{fmpz}, Ref{fmpz_poly}, Int), z, x, 0)
      return ccall((:fmpz_equal, :libflint), Bool,
               (Ref{fmpz}, Ref{fmpz}, Int), z, y, 0)
   else
      return iszero(y)
   end
end

==(x::fmpz, y::fmpz_poly) = y == x

==(x::fmpz_poly, y::Integer) = x == fmpz(y)

==(x::Integer, y::fmpz_poly) = y == x

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(a::fmpz_poly, n::Int)
   n < 0 && throw(DomainError())

   if length(a) <= n
      return a
   end

   z = parent(a)()
   ccall((:fmpz_poly_set_trunc, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, a, n)
   return z
end

function mullow(x::fmpz_poly, y::fmpz_poly, n::Int)
   check_parent(x, y)
   n < 0 && throw(DomainError())

   z = parent(x)()
   ccall((:fmpz_poly_mullow, :libflint), Void,
         (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, x, y, n)
   return z
end

###############################################################################
#
#   Reversal
#
###############################################################################

function reverse(x::fmpz_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpz_poly_reverse, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, x, len)
   return z
end

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpz_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpz_poly_shift_left, :libflint), Void,
      (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, x, len)
   return z
end

function shift_right(x::fmpz_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpz_poly_shift_right, :libflint), Void,
       (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, x, len)
   return z
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   ccall((:fmpz_poly_div, :libflint), Void,
            (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, x, y)
   return z
end

function divrem(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   r = parent(x)()
   ccall((:fmpz_poly_divrem, :libflint), Void,
            (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, r, x, y)
   return z, r
end

function divides(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   flag = Bool(ccall((:fmpz_poly_divides, :libflint), Cint,
           (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, x, y))
   return flag, z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpz_poly, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   ccall((:fmpz_poly_scalar_divexact_fmpz, :libflint), Void,
          (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz}), z, x, y)
   return z
end

function divexact(x::fmpz_poly, y::Int)
   y == 0 && throw(DivideError())
   z = parent(x)()
   ccall((:fmpz_poly_scalar_divexact_si, :libflint), Void,
                        (Ref{fmpz_poly}, Ref{fmpz_poly}, Int), z, x, y)
   return z
end

divexact(x::fmpz_poly, y::Integer) = divexact(x, fmpz(y))

###############################################################################
#
#   Pseudodivision
#
###############################################################################

function pseudorem(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   diff = length(x) - length(y) + 1
   r = parent(x)()
   d = Array{Int}(1)
   ccall((:fmpz_poly_pseudo_rem, :libflint), Void,
     (Ref{fmpz_poly}, Ref{Int}, Ref{fmpz_poly}, Ref{fmpz_poly}), r, d, x, y)
   if (diff > d[1])
      return lead(y)^(diff - d[1])*r
   else
      return r
   end
end

function pseudodivrem(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   diff = length(x) - length(y) + 1
   q = parent(x)()
   r = parent(x)()
   d = Array{Int}(1)
   ccall((:fmpz_poly_pseudo_divrem_divconquer, :libflint), Void,
    (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{Int}, Ref{fmpz_poly}, Ref{fmpz_poly}),
               q, r, d, x, y)
   if (diff > d[1])
      m = lead(y)^(diff - d[1])
      return m*q, m*r
   else
      return q, r
   end
end

###############################################################################
#
#   Content, primitive part, GCD and LCM
#
###############################################################################

function gcd(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpz_poly_gcd, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, x, y)
   return z
end

function content(x::fmpz_poly)
   z = fmpz()
   ccall((:fmpz_poly_content, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_poly}), z, x)
   return z
end

function primpart(x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_primitive_part, :libflint), Void,
         (Ref{fmpz_poly}, Ref{fmpz_poly}), z, x)
   return z
end

###############################################################################
#
#   Evaluation
#
###############################################################################

function evaluate(x::fmpz_poly, y::fmpz)
   z = fmpz()
   ccall((:fmpz_poly_evaluate_fmpz, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz_poly}, Ref{fmpz}), z, x, y)
   return z
end

evaluate(x::fmpz_poly, y::Integer) = evaluate(x, fmpz(y))

###############################################################################
#
#   Composition
#
###############################################################################

function compose(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpz_poly_compose, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, x, y)
   return z
end

###############################################################################
#
#   Derivative
#
###############################################################################

function derivative(x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_derivative, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}), z, x)
   return z
end

###############################################################################
#
#   Resultant
#
###############################################################################

function resultant(x::fmpz_poly, y::fmpz_poly)
   check_parent(x, y)
   z = fmpz()
   ccall((:fmpz_poly_resultant, :libflint), Void,
                (Ref{fmpz}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, x, y)
   return z
end

###############################################################################
#
#   Discriminant
#
###############################################################################

function discriminant(x::fmpz_poly)
   z = fmpz()
   ccall((:fmpz_poly_discriminant, :libflint), Void,
                (Ref{fmpz}, Ref{fmpz_poly}), z, x)
   return z
end

###############################################################################
#
#   RESX
#
###############################################################################

function resx(a::fmpz_poly, b::fmpz_poly)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)
   if lena == 0 || lenb == 0
      return fmpz(), parent(a)(), parent(a)()
   end
   (lena <= 1 && lenb <= 1) && error("Constant polynomials in resx")
   z = fmpz()
   u = parent(a)()
   v = parent(a)()
   c1 = content(a)
   c2 = content(b)
   x = divexact(a, c1)
   y = divexact(b, c2)
   ccall((:fmpz_poly_xgcd_modular, :libflint), Void,
   (Ref{fmpz}, Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}),
            z, u, v, x, y)
   r = z*c1^(lenb - 1)*c2^(lena - 1)
   if lenb > 1
      u *= c1^(lenb - 2)*c2^(lena - 1)
   else
      u *= c2^(lena - 1)
      u = divexact(u, c1)
   end
   if lena > 1
      v *= c1^(lenb - 1)*c2^(lena - 2)
   else
      v *= c1^(lenb - 1)
      v = divexact(v, c2)
   end
   return (r, u, v)
end

###############################################################################
#
#   Signature
#
###############################################################################

doc"""
    signature(f::fmpz_poly)
> Return the signature of the polynomial $f$, i.e. a tuple $(r, s)$ such that
> $r$ is the number of real roots of $f$ and $s$ is half the number of complex
> roots.
"""
function signature(f::fmpz_poly)
   r = Array{Int}(1)
   s = Array{Int}(1)
   ccall((:fmpz_poly_signature, :libflint), Void,
         (Ref{Int}, Ref{Int}, Ref{fmpz_poly}), r, s, f)
   return (r[1], s[1])
end

################################################################################
#
#  Interpolation
#
################################################################################

function interpolate(R::FmpzPolyRing, x::Array{fmpz, 1},
                                      y::Array{fmpz, 1})
  z = R()

  ax = Array{Int}(length(x))
  ay = Array{Int}(length(y))

  t = fmpz()

  for i in 1:length(x)
    ax[i] = x[i].d
    ay[i] = y[i].d
  end

  ccall((:fmpz_poly_interpolate_fmpz_vec, :libflint), Void,
          (Ref{fmpz_poly}, Ref{Int}, Ref{Int}, Int),
          z, ax, ay, length(x))
  return z
end

################################################################################
#
#  Factorization
#
################################################################################

doc"""
    factor(x::fmpz_poly)
> Returns the factorization of $x$.
"""
function factor(x::fmpz_poly)
  fac, z = _factor(x)
  ffac = factor(z)

  for (p, e) in ffac
    fac[parent(x)(p)] = e
  end

  return Fac(parent(x)(unit(ffac)), fac)
end

function _factor(x::fmpz_poly)
  fac = fmpz_poly_factor()
  ccall((:fmpz_poly_factor, :libflint), Void,
              (Ref{fmpz_poly_factor}, Ref{fmpz_poly}), fac, x)
  res = Dict{fmpz_poly,Int}()
  z = fmpz()
  ccall((:fmpz_poly_factor_get_fmpz, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz_poly_factor}), z, fac)
  for i in 1:fac.num
    f = parent(x)()
    ccall((:fmpz_poly_factor_get_fmpz_poly, :libflint), Void,
            (Ref{fmpz_poly}, Ref{fmpz_poly_factor}, Int), f, fac, i - 1)
    e = unsafe_load(fac.exp, i)
    res[f] = e
  end
  return res, z
end

###############################################################################
#
#   Special polynomials
#
###############################################################################

function chebyshev_t(n::Int, x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_chebyshev_t, :libflint), Void,
                                                  (Ref{fmpz_poly}, Int), z, n)
   return isgen(x) ? z : compose(z, x)
end

function chebyshev_u(n::Int, x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_chebyshev_u, :libflint), Void,
                                                  (Ref{fmpz_poly}, Int), z, n)
   return isgen(x) ? z : compose(z, x)
end

doc"""
    cyclotomic(n::Int, x::fmpz_poly)
> Return the $n$th cyclotomic polynomial, defined as
> $$\Phi_n(x) = \prod_{\omega} (x-\omega),$$ where $\omega$ runs over all the
> $n$th primitive roots of unity.
"""
function cyclotomic(n::Int, x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_cyclotomic, :libflint), Void,
                                                  (Ref{fmpz_poly}, Int), z, n)
   return isgen(x) ? z : compose(z, x)
end

doc"""
    swinnerton_dyer(n::Int, x::fmpz_poly)
> Return the Swinnerton-Dyer polynomial $S_n$, defined as the integer
> polynomial
> $$S_n = \prod (x \pm \sqrt{2} \pm \sqrt{3} \pm \sqrt{5} \pm \ldots \pm \sqrt{p_n})$$
> where $p_n$ denotes the $n$-th prime number and all combinations of signs are
> taken. This polynomial has degree $2^n$ and is irreducible over the integers
> (it is the minimal polynomial of $\sqrt{2} + \ldots + \sqrt{p_n}$).
"""
function swinnerton_dyer(n::Int, x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_swinnerton_dyer, :libflint), Void,
                                                  (Ref{fmpz_poly}, Int), z, n)
   return isgen(x) ? z : compose(z, x)
end

doc"""
    cos_minpoly(n::Int, x::fmpz_poly)
> Return the minimal polynomial of $2 \cos(2 \pi / n)$. For suitable choice of
> $n$, this gives the minimal polynomial of $2 \cos(a \pi)$ or $2 \sin(a \pi)$ for any
> rational $a$.
"""
function cos_minpoly(n::Int, x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_cos_minpoly, :libflint), Void,
                                                  (Ref{fmpz_poly}, Int), z, n)
   return isgen(x) ? z : compose(z, x)
end

doc"""
    theta_qexp(e::Int, n::Int, x::fmpz_poly)
> Return the $q$-expansion to length $n$ of the Jacobi theta function raised to
> the power $r$, i.e. $\vartheta(q)^r$ where
> $\vartheta(q) = 1 + \sum_{k=1}^{\infty} q^{k^2}$.
"""
function theta_qexp(e::Int, n::Int, x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_theta_qexp, :libflint), Void,
                                          (Ref{fmpz_poly}, Int, Int), z, e, n)
   return isgen(x) ? z : compose(z, x)
end

doc"""
    eta_qexp(e::Int, n::Int, x::fmpz_poly)
> Return the $q$-expansion to length $n$ of the Dedekind eta function (without
> the leading factor $q^{1/24}$) raised to the power $r$, i.e.
> $(q^{-1/24} \eta(q))^r = \prod_{k=1}^{\infty} (1 - q^k)^r$.
> In particular, $r = -1$ gives the generating function of the partition
> function $p(k)$, and $r = 24$ gives, after multiplication by $q$, the modular
> discriminant $\Delta(q)$ which generates the Ramanujan tau function
> $\tau(k)$.
"""
function eta_qexp(e::Int, n::Int, x::fmpz_poly)
   z = parent(x)()
   ccall((:fmpz_poly_eta_qexp, :libflint), Void,
                                          (Ref{fmpz_poly}, Int, Int), z, e, n)
   return isgen(x) ? z : compose(z, x)
end

###############################################################################
#
#   Speedups for polynomials over fmpz_polys
#
###############################################################################

function *(a::Generic.Poly{fmpz_poly}, b::Generic.Poly{fmpz_poly})
   check_parent(a, b)
   if min(length(a), length(b)) < 40
      return mul_classical(a, b)
   else
      return mul_ks(a, b)
   end
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(z::fmpz_poly)
   ccall((:fmpz_poly_zero, :libflint), Void,
                    (Ref{fmpz_poly},), z)
   return z
end

function fit!(z::fmpz_poly, n::Int)
   ccall((:fmpz_poly_fit_length, :libflint), Void,
                    (Ref{fmpz_poly}, Int), z, n)
   return nothing
end

function setcoeff!(z::fmpz_poly, n::Int, x::fmpz)
   ccall((:fmpz_poly_set_coeff_fmpz, :libflint), Void,
                    (Ref{fmpz_poly}, Int, Ref{fmpz}), z, n, x)
   return z
end

function mul!(z::fmpz_poly, x::fmpz_poly, y::fmpz_poly)
   ccall((:fmpz_poly_mul, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, x, y)
   return z
end

function addeq!(z::fmpz_poly, x::fmpz_poly)
   ccall((:fmpz_poly_add, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, z, x)
   return z
end

function add!(z::fmpz_poly, x::fmpz_poly, y::fmpz_poly)
   ccall((:fmpz_poly_add, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}), z, x, y)
   return z
end

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{fmpz_poly}, ::Type{T}) where {T <: Integer} = fmpz_poly

promote_rule(::Type{fmpz_poly}, ::Type{fmpz}) = fmpz_poly

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (a::FmpzPolyRing)()
   z = fmpz_poly()
   z.parent = a
   return z
end

function (a::FmpzPolyRing)(b::Int)
   z = fmpz_poly(b)
   z.parent = a
   return z
end

function (a::FmpzPolyRing)(b::Integer)
   z = fmpz_poly(fmpz(b))
   z.parent = a
   return z
end

function (a::FmpzPolyRing)(b::fmpz)
   z = fmpz_poly(b)
   z.parent = a
   return z
end

function (a::FmpzPolyRing)(b::Array{fmpz, 1})
   z = fmpz_poly(b)
   z.parent = a
   return z
end

(a::FmpzPolyRing)(b::Array{T, 1}) where {T <: Integer} = a(map(fmpz, b))

(a::FmpzPolyRing)(b::fmpz_poly) = b

###############################################################################
#
#   PolynomialRing constructor
#
###############################################################################

function PolynomialRing(R::FlintIntegerRing, s::AbstractString; cached = true)
   S = Symbol(s)

   parent_obj = FmpzPolyRing(S, cached)

   return parent_obj, parent_obj([fmpz(0), fmpz(1)])
end


text:
#=

  Set covering and set partition in Julia ConstraintSolver.jl

  Example from Lundgren, Ronnqvist, Varbrand "Optimeringslora", page 408.
  [This is a Swedish book about Operational Research.]

  We want to minimize the cost of the alternatives which covers all the
  objects, i.e. all objects must be choosen. The requirement is than an object
  may be selected _exactly_ once.

  Alternative        Cost        Object
  1                  19           1,6
  2                  16           2,6,8
  3                  18           1,4,7
  4                  13           2,3,5
  5                  15           2,5
  6                  19           2,3
  7                  15           2,3,4
  8                  17           4,5,8
  9                  16           3,6,8
  10                 15           1,6,7

  The problem has a unique solution of z = 49 where alternatives
  3, 5, and 9 is selected.

  If we, however, allow that an object is selected more than one time,
  then the solution is z = 45 (i.e. less cost than the first problem),
  and the alternatives 4, 8, and 10 is selected, where object 5 is
  selected twice (alt. 4 and 8). It's an unique solution as well.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function set_covering4(problem,type="",print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 100, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    println("Type:$type")
    costs = problem[:costs]
    alternatives = problem[:alternatives]

    num_alternatives, num_objects = size(alternatives)

    # Which alternative to select
    @variable(model, x[1:num_alternatives], Bin)
    @variable(model, 0 <= min_val <= num_alternatives*maximum(costs), Int)

    for j in 1:num_objects
        s = @variable(model, [1:1], CS.Integers(0:num_alternatives*sum(alternatives[:,j])))
        @constraint(model, s[1] == sum(x.*alternatives[:,j]))
        if type == "set_partition"
            # all objects must be covered _exactly_ once
            # (set partition)
            @constraint(model, s[1] == 1)
        else
            # variant: all objects must be covered _at least_ once
            # (set covering)
            @constraint(model, s[1] >= 1)
        end
    end

    scalar_product(model,costs,x,min_val)

    @objective(model,Min,min_val)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                min_val_val = convert.(Integer,JuMP.value.(min_val; result=sol))
                println("min_val:$min_val_val")
                println("x:$x_val")
                println("assignments:$([i for i in 1:num_alternatives if x_val[i] == 1])")
            end
        end
    else
        println("status:$status")
    end

    return status
end

problem = Dict(
    :costs => [9, 16, 18, 13, 15, 19, 15, 17, 16, 15],
    :alternatives => resize_matrix([
        [1,0,0,0,0,1,0,0],   # alternative 1
        [0,1,0,0,0,1,0,1],   # alternative 2
        [1,0,0,1,0,0,1,0],   # alternative 3
        [0,1,1,0,1,0,0,0],   # alternative 4
        [0,1,0,0,1,0,0,0],   # alternative 5
        [0,1,1,0,0,0,0,0],   # alternative 6
        [0,1,1,1,0,0,0,0],   # alternative 7
        [0,0,0,1,1,0,0,1],   # alternative 8
        [0,0,1,0,0,1,0,1],   # alternative 9
        [1,0,0,0,0,1,1,0]]  # alternative 10
                                   )
)

@time set_covering4(problem,"set_partition",true,true)
println()
@time set_covering4(problem,"set_covering",true,true)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module Helper

"""# UnExported

	@console 	expresion 	[console_inputs...]

execute the expresion `expresion` after putting `console_inputs` into the standar input buffer.

for functions that meant to run iteractivelly while require user intput, this macro provid a tool to
allow for noninteractive testing such functions by providing the input in advance by `console_inputs`.
"""
macro console(expression, console_inputs...)
	quote
		bffr = readavailable(stdin.buffer) # empty input stream to ensure later only the `console_inputs` is in `stdin` buffer.
		bffr == UInt8[] || @warn "buffer not empty, see the pervious process to 'stdin'"  buffer = String(bffr)
		if 0 == length($console_inputs)
			write(stdin.buffer, "\n")   					# [1of2] empty `console_inputs` simulate return or enter .
		else
			for input in $console_inputs
				if typeof(input) == String && input != "" 	# [2of2]empty `console_inputs` simulate return or enter .
					for npt in string.(split(input))
						write(stdin.buffer, npt, "\n")
					end
				else
					write(stdin.buffer, string(input), "\n")
				end #if
			end # for
		end #IF
		$expression		#:($(esc(expression)))
	end |> esc
end

"""# UnExported

	exec_console_unattended(Fn::Union{Function, Type}, console_inputs::Vector = [];
																	Fn_ARGs::Vector=[])

execute the expresion `Fn(Fn_ARGs...)` after putting `console_inputs` into the standar input buffer.

for functions that meant to run iteractivelly while require user intput, this function provid a tool to
allow for noninteractive testing such functions by providing the input in advance by `console_inputs`.
"""
function exec_console_unattended(Fn::Union{Function, Type}, console_inputs::Vector = []; Fn_ARGs::Vector=[])
	bffr = readavailable(stdin.buffer) # empty input stream to ensure later only the `console_inputs` is in `stdin` buffer.
	bffr == UInt8[] || @warn "buffer not empty, see the pervious process to 'stdin'"  buffer = String(bffr)

	if 0 == length(console_inputs)
		write(stdin.buffer, "\n")   # empty `console_inputs` simulate return or enter .

	else
		for input in  string.(console_inputs)
			write(stdin.buffer, input, "\n")
		end # for

	end #IF

	return Fn(Fn_ARGs...)		# call and return the value
end

"""# UnExported

	exec_console_unattended(Fn::Union{Function, Type}, console_inputs...;
																	Fn_ARGs::Vector=[])

execute the expresion `Fn(Fn_ARGs...)` after putting `console_inputs` into the standar input buffer.

for functions that meant to run iteractivelly while require user intput, this function provid a tool to
allow for noninteractive testing such functions by providing the input in advance by `console_inputs`.
"""
exec_console_unattended(Fn::Union{Function, Type}, console_inputs...; Fn_ARGs::Vector=[]) = exec_console_unattended(Fn, [console_inputs...]; Fn_ARGs=Fn_ARGs)

"""# UnExported

	exec_console_unattended(Fn::Union{Function, Type}, console_inputs::String;
																	Fn_ARGs::Vector=[])

execute the expresion `Fn(Fn_ARGs...)` after putting `console_inputs` into the standar input buffer.

for functions that meant to run iteractivelly while require user intput, this function provid a tool to
allow for noninteractive testing such functions by providing the input in advance by `console_inputs`.
"""
exec_console_unattended(Fn::Union{Function, Type}, console_inputs::String; Fn_ARGs::Vector=[]) = exec_console_unattended(Fn, split(console_inputs); Fn_ARGs=Fn_ARGs)

end
const H = Helper


text:
module PolicyCompliance

include(joinpath("Units.jl"))

end


text:
#TODO: add @inline, test with np >= 2, test with CuArrays
using Base.Broadcast
import Base.Broadcast: BroadcastStyle, Broadcasted
import Base: tail

struct MPIArrayStyle{Style <: BroadcastStyle} <: Broadcast.AbstractArrayStyle{Any} end
MPIArrayStyle(::S) where {S} = MPIArrayStyle{S}()
MPIArrayStyle(::S, ::Val{N}) where {S,N} = MPIArrayStyle(S(Val(N)))
MPIArrayStyle(::Val{N}) where {N} = MPIArrayStyle{Broadcast.DefaultArrayStyle{N}}()

BroadcastStyle(::Type{<:MPIArray{T, N, A}}) where {T,N,A} = MPIArrayStyle(BroadcastStyle(A{T,N}), Val(N))

function BroadcastStyle(::MPIArrayStyle{AStyle}, ::MPIArrayStyle{BStyle}) where {AStyle, BStyle}
    MPIArrayStyle(BroadcastStyle(AStyle, BStyle))()
end

BroadcastStyle(::Type{<:SubArray{<:Any, <:Any, <:T}}) where T <: MPIArray = BroadcastStyle(T)

function Broadcast.broadcasted(::MPIArrayStyle{Style}, f, args...) where Style
    inner = Broadcast.broadcasted(Style(), f, args...)
    if inner isa Broadcasted
        return Broadcasted{MPIArrayStyle{Style}}(inner.f, inner.args, inner.axes)
    else
        return inner
    end
end

# TODO: For Transpose, Adjoint, SubArray, etc.
const MPIDestArray = MPIArray

function Base.similar(bc::Broadcasted{<:MPIArrayStyle{Style}}, ::Type{ElType}) where {Style, ElType}
    MPIArray(map(length, axes(bc))...) do I
        bcc = Broadcasted{Style}(identity, (), map(length, I))
        similar(bcc, ElType)
    end
end

function Base.copyto!(dest::MPIDestArray, bc::Broadcasted)
    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
    dbc = bcdistribute(bc)

    lidcs = dest.partitioning[Rank()+1]
    lbc = bclocal(dbc, lidcs)
    Base.copyto!(dest.localarray, lbc)
    return dest
end

function Base.copy(bc::Broadcasted{<:MPIArrayStyle})
    dbc = bcdistribute(bc)
    MPIArray(map(length, axes(bc))...) do I
        lbc = Broadcast.instantiate(bclocal(dbc, I))
        copy(lbc)
    end
end

bcdistribute(bc::Broadcasted{Style}) where Style = Style <: BroadcastStyle ? Broadcasted{MPIArrayStyle{Style}}(bc.f, bcdistribute_args(bc.args), bc.axes) : bc
bcdistribute(bc::Broadcasted{Style}) where Style <: MPIArrayStyle  = Broadcasted{Style}(bc.f, bcdistribute_args(bc.args), bc.axes)

bcdistribute(x::T) where T = _bcdistribute(BroadcastStyle(T), x)
_bcdistribute(::MPIArrayStyle, x) = x
_bcdistribute(::Broadcast.AbstractArrayStyle{0}, x) = x
_bcdistribute(::Broadcast.AbstractArrayStyle, x) = x # error("not implemented") # distribute(x) TODO: define distribute
_bcdistribute(::Broadcast.AbstractArrayStyle, x::AbstractRange) = x
_bcdistribute(::Any, x) = x

bcdistribute_args(args::Tuple) = (bcdistribute(args[1]), bcdistribute_args(tail(args))...)
bcdistribute_args(args::Tuple{Any}) = (bcdistribute(args[1]),)
bcdistribute_args(args::Tuple{}) = ()

_bcview(::Tuple{}, ::Tuple{}) = ()
_bcview(::Tuple{}, view::Tuple) = ()
_bcview(shape::Tuple, ::Tuple{}) = (shape[1], _bcview(tail(shape),())...)
function _bcview(shape::Tuple, view::Tuple)
    return (_bcview1(shape[1], view[1]), _bcview(tail(shape), tail(view))...)
end

# _bcview1 handles the logic for a single dimension
function _bcview1(a, b)
    if a==1 || a == 1:1
        return 1:1
    elseif first(a) <= first(b) <= last(a) &&
            first(a) <= last(b) <= last(b)
        return b
    else
        throw(DimensionMismatch("broadcast view could not be constructed"))
    end
end

@inline bclocal(bc::Broadcasted{MPIArrayStyle{Style}}, idxs) where Style = Broadcasted{Style}(bc.f, bclocal_args(_bcview(axes(bc), idxs), bc.args))

# bclocal(bc::Broadcasted{Nothing}, idxs) = bc
bclocal(x::T, idxs) where T = _bclocal(BroadcastStyle(T), x, idxs)
function _bclocal(::MPIArrayStyle, x, idxs)
    bcidxs = _bcview(axes(x), idxs)
    @assert all(map( (y, z) -> (y  == z) || length(z) == 1 , bcidxs, x.partitioning[Rank() + 1]))
    x.localarray
    # makelocal(x, bcidxs...) # TODO: makelocal
end

_bclocal(::Broadcast.AbstractArrayStyle{0}, x, idxs) = x
function _bclocal(::Broadcast.AbstractArrayStyle, x::AbstractRange, idxs)
    @assert length(idxs) == 1
    x[idxs[1]]
end
function _bclocal(::Broadcast.Style{Tuple}, x, idxs)
    @assert length(idxs) == 1
    tuple((e for (i,e) in enumerate(x) if i in idxs[1])...)
end
_bclocal(::Broadcast.AbstractArrayStyle, x, idxs) = x
_bclocal(::Any, x, idxs) = error("don't know how to localize $x with $idxs")

@inline bclocal_args(idxs, args::Tuple) = (bclocal(args[1], idxs), bclocal_args(idxs, tail(args))...)
bclocal_args(idxs, args::Tuple{Any}) = (bclocal(args[1], idxs),)
bclocal_args(idxs, args::Tuple{}) = ()


text:
using ShowEach
using Base.Test

macro test_output(expr, output)
    quote
        origSTDOUT = STDOUT
        (rd, wr) = redirect_stdout()
        ShowEach.state.io = wr
        $(esc(expr))
        print(" ") # work around blocking of `readavailable` for empty pipe
        close(wr)
        expr_output = chop(String(readavailable(rd)))
        close(rd)
        redirect_stdout(origSTDOUT)
        ShowEach.state.io = origSTDOUT
        @test expr_output == $output
    end
end

# Some definitions for testing
type A{T} x::T end
a = [1,2,3]
x = 1
y = 2
z = 3
f(x, y) = x * y
g(z) = 2z
o = A(3)

set_loop_limit(20)

@test_output (@showeach 1) ""
@test_output (@showeach x) "x = 1\n"
@test_output (@showeach a += x) "x = 1\n"
@test_output (@showeach x + 1) """
┌─ x = 1
x + 1 = 2
"""
@test_output (@showeach a + 1) """
┌─ a = [2,3,4]
a + 1 = [3,4,5]
"""
@test_output (@showeach f(x, y)) """
┌─ x = 1
├─ y = 2
f(x,y) = 2
"""
@test_output (@showeach z*f(x,y)) """
┌─ z = 3
│  ┌─ x = 1
│  ├─ y = 2
├─ f(x,y) = 2
z * f(x,y) = 6
"""
@test_output (@showeach f(g(x), a[y])) """
   ┌─ x = 1
┌─ g(x) = 2
│  ┌─ a = [2,3,4]
│  ├─ y = 2
├─ a[y] = 3
f(g(x),a[y]) = 6
"""
@test_output (@showeach f.(a, [x,y,z])) """
┌─ a = [2,3,4]
├─ x = 1
├─ y = 2
├─ z = 3
f.(a,[x,y,z]) = [2,6,12]
"""
@test_output (@showeach o) """
o = A{Int64}(3)
"""
@test_output (@showeach g(o.x)) """
┌─ o.x = 3
g(o.x) = 6
"""
@test_output (@showeach g(A{Float64}(1).x)) """
┌─ A{Float64}(1).x = 1.0
g(A{Float64}(1).x) = 2.0
"""

myrand(x) = [0.159358,0.35803,0.791544]
@test_output (@showeach g.(myrand(3)[2])) """
   ┌─ myrand(3) = [0.159358,0.35803,0.791544]
┌─ (myrand(3))[2] = 0.35803
g.((myrand(3))[2]) = 0.71606
"""
@test_output (@showeach Dict(:a => 2, :b =>  3)[:a]) """
┌─ Dict(:a => 2,:b => 3) = Dict{Symbol,Int64} with 2 entries:
│                            :a => 2
│                            :b => 3
(Dict(:a => 2,:b => 3))[:a] = 2
"""
@test_output (@showeach sum(i for i = x:y if iseven(i))) """
┌─ x = 1
├─ y = 2
│  ┌─ i = 1
├─ iseven(i) = false
│  ┌─ i = 2
├─ iseven(i) = true
├─ i = 2
sum((i for i = x:y if iseven(i))) = 2
"""
@test_output (@showeach sum(i for i = x:y for j = i:z if iseven(i))) """
┌─ x = 1
├─ y = 2
├─ i = 1
├─ z = 3
│  ┌─ i = 1
├─ iseven(i) = false
│  ┌─ i = 1
├─ iseven(i) = false
│  ┌─ i = 1
├─ iseven(i) = false
├─ i = 2
├─ z = 3
│  ┌─ i = 2
├─ iseven(i) = true
│  ┌─ i = 2
├─ iseven(i) = true
├─ i = 2
├─ i = 2
sum((i for i = x:y for j = i:z if iseven(i))) = 4
"""
@test_output (@showeach A{Float64}) """
"""
@test_output (@showeach A{Float64}(1)) """
A{Float64}(1) = A{Float64}(1.0)
"""
@test_output (@showeach [i*j for i in x:y, j in eachindex(a)]) """
┌─ x = 1
├─ y = 2
│  ┌─ a = [2,3,4]
├─ eachindex(a) = Base.OneTo(3)
│  ┌─ i = 1
│  ├─ j = 1
├─ i * j = 1
│  ┌─ i = 2
│  ├─ j = 1
├─ i * j = 2
│  ┌─ i = 1
│  ├─ j = 2
├─ i * j = 2
│  ┌─ i = 2
│  ├─ j = 2
├─ i * j = 4
│  ┌─ i = 1
│  ├─ j = 3
├─ i * j = 3
│  ┌─ i = 2
│  ├─ j = 3
├─ i * j = 6
[i * j for i = x:y, j = eachindex(a)] = 2×3 Array{Int64,2}:
                                            1  2  3
                                            2  4  6
"""
@test_output (@showeach Float64[i*j for i in x:y, j in eachindex(a)]) """
┌─ x = 1
├─ y = 2
│  ┌─ a = [2,3,4]
├─ eachindex(a) = Base.OneTo(3)
│  ┌─ i = 1
│  ├─ j = 1
├─ i * j = 1
│  ┌─ i = 2
│  ├─ j = 1
├─ i * j = 2
│  ┌─ i = 1
│  ├─ j = 2
├─ i * j = 2
│  ┌─ i = 2
│  ├─ j = 2
├─ i * j = 4
│  ┌─ i = 1
│  ├─ j = 3
├─ i * j = 3
│  ┌─ i = 2
│  ├─ j = 3
├─ i * j = 6
Float64[i * j for i = x:y, j = eachindex(a)] = 2×3 Array{Float64,2}:
                                                   1.0  2.0  3.0
                                                   2.0  4.0  6.0
"""
@test_output (@showeach a[1] < x < b) """
┌─ x = 1
│  ┌─ a = [2,3,4]
├─ a[1] = 2
a[1] < x < b = false
"""
@test_output (@showeach f(x) = 2x) """
"""
@test_output f(2) """
in function f(x):
   ┌─ x = 2
   2x = 4
"""
@test_output (@showeach function h(x)
    a = 2x
    println("foo")
    a
end) """
"""
@test_output (h(x)) """
in function h(x):
   ┌─ x = 1
   2x = 2
foo
   println("foo") = nothing
   a = 2
"""
@test_output (@showeach f(x...) = x) """
"""
@test_output (@showeach f(a...)) """
┌─ a = [2,3,4]
├─ in function f(x...):
│     x = (2,3,4)
f(a...) = (2,3,4)
"""
@test_output (@showeach if x < 2 a else b end) """
┌─ x = 1
x < 2 = true
a = [2,3,4]
"""
# TODO: let blocks
@test_output (@showeach let x = 1
    y += 1
    x + y
end) """
"""
y = 3
@test_output (@showeach for i in x:y, j = a
            i*j
        end) """
┌─ x = 1
├─ y = 3
├─ a = [2,3,4]
│  ┌─ i = 1
│  ├─ j = 2
├─ i * j = 2
│  ┌─ i = 1
│  ├─ j = 3
├─ i * j = 3
│  ┌─ i = 1
│  ├─ j = 4
├─ i * j = 4
├─ a = [2,3,4]
│  ┌─ i = 2
│  ├─ j = 2
├─ i * j = 4
│  ┌─ i = 2
│  ├─ j = 3
├─ i * j = 6
│  ┌─ i = 2
│  ├─ j = 4
├─ i * j = 8
├─ a = [2,3,4]
│  ┌─ i = 3
│  ├─ j = 2
├─ i * j = 6
│  ┌─ i = 3
│  ├─ j = 3
├─ i * j = 9
│  ┌─ i = 3
│  ├─ j = 4
├─ i * j = 12
for ... end = nothing
"""
y = 2
@test_output (@showeach x < a[3] || error("Foo")) """
   ┌─ x = 1
   │  ┌─ a = [2,3,4]
   ├─ a[3] = 4
┌─ x < a[3] = true
x < a[3] || error("Foo") = true
"""
@test_output (@showeach x < a[3] && (x += 1)) """
   ┌─ x = 1
   │  ┌─ a = [2,3,4]
   ├─ a[3] = 4
┌─ x < a[3] = true
x < a[3] && (x += 1) = 2
"""
x = 1

@test_output (@showeach while x < z
    x += 1
end) """
   ┌─ x = 1
   ├─ z = 3
┌─ x < z = true
│  ┌─ x = 2
│  ├─ z = 3
├─ x < z = true
│  ┌─ x = 3
│  ├─ z = 3
├─ x < z = false
while ... end = nothing
"""
x = 1

@test_output (@showeach g(x)) """
┌─ x = 1
g(x) = 2
"""
@test_output (@showeach f(x,y)) """
┌─ x = 1
├─ y = 2
f(x,y) = 2
"""
@test_output (@showeach f(g(y), g(x))) """
   ┌─ y = 2
┌─ g(y) = 4
│  ┌─ x = 1
├─ g(x) = 2
f(g(y),g(x)) = 8
"""
@test_output (@showeach f(f(x), f(y))) """
   ┌─ x = 1
   ├─ in function f(x):
   │     ┌─ x = 1
   │     2x = 2
┌─ f(x) = 2
│  ┌─ y = 2
│  ├─ in function f(x):
│  │     ┌─ x = 2
│  │     2x = 4
├─ f(y) = 4
f(f(x),f(y)) = 8
"""
@test_output (@showeach function p(x, y)
           println("foo")
           ret = x*y
           println("bar")
           ret
       end) """
"""
a = [1,2,3]
@test_output (@showeach begin
       println("enter block")
       a*p(x,y)
       println("in block")
       z*p(x,y)
       println("end block")
   end) """
enter block
println("enter block") = nothing
┌─ a = [1,2,3]
│  ┌─ x = 1
│  ├─ y = 2
│  ├─ in function p(x,y):
foo
│  │     println("foo") = nothing
│  │     ┌─ x = 1
│  │     ├─ y = 2
│  │     x * y = 2
bar
│  │     println("bar") = nothing
│  │     ret = 2
├─ p(x,y) = 2
a * p(x,y) = [2,4,6]
in block
println("in block") = nothing
┌─ z = 3
│  ┌─ x = 1
│  ├─ y = 2
│  ├─ in function p(x,y):
foo
│  │     println("foo") = nothing
│  │     ┌─ x = 1
│  │     ├─ y = 2
│  │     x * y = 2
bar
│  │     println("bar") = nothing
│  │     ret = 2
├─ p(x,y) = 2
z * p(x,y) = 6
end block
println("end block") = nothing
"""

gather_stdout(true)
@test_output (@showeach begin
       println("enter block")
       a*p(x,y)
       println("in block")
       z*p(x,y)
       println("end block")
   end) """
println("enter block") = nothing
┌─ a = [1,2,3]
│  ┌─ x = 1
│  ├─ y = 2
│  ├─ in function p(x,y):
│  │     println("foo") = nothing
│  │     ┌─ x = 1
│  │     ├─ y = 2
│  │     x * y = 2
│  │     println("bar") = nothing
│  │     ret = 2
├─ p(x,y) = 2
a * p(x,y) = [2,4,6]
println("in block") = nothing
┌─ z = 3
│  ┌─ x = 1
│  ├─ y = 2
│  ├─ in function p(x,y):
│  │     println("foo") = nothing
│  │     ┌─ x = 1
│  │     ├─ y = 2
│  │     x * y = 2
│  │     println("bar") = nothing
│  │     ret = 2
├─ p(x,y) = 2
z * p(x,y) = 6
println("end block") = nothing
gathered STDOUT:
enter block
foo
bar
in block
foo
bar
end block
"""

set_loop_limit(6)
@test_output (@showeach p(x) = 2x) """
"""
@test_output (@showeach map(p, 1:100)) """
┌─ p = p (generic function with 2 methods)
├─ in function p(x):
│     ┌─ x = 1
│     2x = 2
├─ in function p(x):
│     ┌─ x = 2
│     2x = 4
├─ in function p(x):
│     ┌─ x = 3
│     2x = 6
├─ in function p(x):
│     ┌─ x = 4
│     2x = 8
├─ in function p(x):
│     ┌─ x = 5
│     2x = 10
├─  ⋮
map(p,1:100) = [2,4,6  …  196,198,200]
"""
@test_output (@showeach sum(i for i = 1:10 if isodd(i))) """
   ┌─ i = 1
┌─ isodd(i) = true
│  ┌─ i = 2
├─ isodd(i) = false
│  ┌─ i = 3
├─ isodd(i) = true
├─ i = 1
│  ┌─ i = 4
├─ isodd(i) = false
│  ┌─ i = 5
├─ isodd(i) = true
├─ i = 3
├─  ⋮
sum((i for i = 1:10 if isodd(i))) = 25
"""