text:
# Test basic extensions
@testset "Basics" begin
    # initialize the model, references, and other info
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 10])
    con1 = ScalarConstraint(zero(AffExpr), MOI.EqualTo(0.0))
    con2 = BoundedScalarConstraint(zero(AffExpr), MOI.EqualTo(0.0),
                                  ParameterBounds((t => IntervalSet(0, 1),)),
                                  ParameterBounds())
    object1 = ConstraintData(con1, Int[], "c1", MeasureIndex[], false)
    object2 = ConstraintData(con2, [1], "c2", MeasureIndex[], false)
    idx1 = ConstraintIndex(1)
    idx2 = ConstraintIndex(2)
    cref1 = InfOptConstraintRef(m, idx1)
    cref2 = InfOptConstraintRef(m, idx2)
    bad_cref = InfOptConstraintRef(m, ConstraintIndex(-1))
    # owner_model
    @testset "JuMP.owner_model" begin
      @test owner_model(cref1) === m
      @test owner_model(cref2) === m
    end
    # index
    @testset "JuMP.index" begin
      @test JuMP.index(cref1) == ConstraintIndex(1)
      @test JuMP.index(cref2) == ConstraintIndex(2)
    end
    # test Base.:(==) of constraint references
    @testset "Base.:(==) References" begin
        @test !(cref1 == cref2)
        @test cref1 == InfOptConstraintRef(m, idx1)
        @test cref1 != InfOptConstraintRef(m, idx2)
        @test cref1 != InfOptConstraintRef(InfiniteModel(), idx1)
    end
    # test broadcastable
    @testset "Base.broadcastable Reference" begin
        @test isa(Base.broadcastable(cref1), Base.RefValue)
    end
    # test constraint type
    @testset "JuMP.constraint_type" begin
        @test constraint_type(m) == InfOptConstraintRef
    end
    # test shape
    @testset "JuMP.shape" begin
        @test JuMP.shape(con2) == ScalarShape()
    end
    # test jump_function
    @testset "JuMP.jump_function" begin
        @test jump_function(con2) == zero(AffExpr)
    end
    # test moi_set
    @testset "JuMP.moi_set" begin
        @test moi_set(con2) == MOI.EqualTo(0.0)
    end
    # test parameter_bounds
    @testset "parameter_bounds" begin
        @test parameter_bounds(con2) == ParameterBounds((t => IntervalSet(0, 1),))
    end
    # test original_parameter_bounds
    @testset "original_parameter_bounds" begin
        @test original_parameter_bounds(con2) == ParameterBounds()
    end
    # _add_data_object
    @testset "_add_data_object" begin
        @test InfiniteOpt._add_data_object(m, object1) == idx1
        @test InfiniteOpt._add_data_object(m, object2) == idx2
    end
    # _data_dictionary
    @testset "_data_dictionary" begin
        @test InfiniteOpt._data_dictionary(cref1) === m.constraints
    end
    # is_valid
    @testset "JuMP.is_valid" begin
      @test is_valid(m, cref1)
      @test is_valid(m, cref2)
    end
    # _data_object
    @testset "_data_object" begin
        @test InfiniteOpt._data_object(cref1) === object1
        @test InfiniteOpt._data_object(cref2) === object2
        @test_throws ErrorException InfiniteOpt._data_object(bad_cref)
    end
    # _core_constraint_object
    @testset "_core_constraint_object" begin
        @test InfiniteOpt._core_constraint_object(cref1) === con1
        @test InfiniteOpt._core_constraint_object(cref2) === con2
    end
    # _set_core_constraint_object
    @testset "_set_core_constraint_object" begin
        @test InfiniteOpt._set_core_constraint_object(cref1, con2) isa Nothing
        @test InfiniteOpt._set_core_constraint_object(cref1, con1) isa Nothing
    end
    # _object_numbers
    @testset "_object_numbers" begin
        @test InfiniteOpt._object_numbers(cref1) == []
        @test InfiniteOpt._object_numbers(cref2) == [1]
    end
    # _measure_dependencies
    @testset "_measure_dependencies" begin
        @test InfiniteOpt._measure_dependencies(cref1) == MeasureIndex[]
    end
    @testset "_is_info_constraint" begin
        @test !InfiniteOpt._is_info_constraint(cref2)
    end
    # constraint_object
    @testset "JuMP.constraint_object" begin
      @test constraint_object(cref1) == con1
      @test constraint_object(cref2) == con2
    end
    # JuMP.name
    @testset "JuMP.name" begin
        @test name(cref1) == "c1"
        @test name(cref2) == "c2"
        @test name(bad_cref) == ""
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        @test isa(set_name(cref1, "new"), Nothing)
        @test_throws ErrorException set_name(bad_cref, "")
        @test name(cref1) == "new"
    end
    # test has_parameter_bounds
    @testset "has_parameter_bounds" begin
        @test has_parameter_bounds(cref2)
        @test !has_parameter_bounds(cref1)
    end
    # test parameter_bounds
    @testset "parameter_bounds" begin
        @test_throws ErrorException parameter_bounds(cref1)
        @test parameter_bounds(cref2) == ParameterBounds((t => IntervalSet(0, 1),))
    end
    # _make_constraint_ref
    @testset "_make_constraint_ref" begin
        @test InfiniteOpt._make_constraint_ref(m, idx1) == cref1
        @test InfiniteOpt._make_constraint_ref(m, idx2) == cref2
    end
    # constraint_by_name
    @testset "JuMP.constraint_by_name" begin
        # test normal
        @test constraint_by_name(m, "new") == cref1
        @test isa(constraint_by_name(m, "test2"), Nothing)
        # prepare constraint with duplicate name
        @test isa(set_name(cref2, "new"), Nothing)
        # test for duplciate name error
        @test_throws ErrorException constraint_by_name(m, "new")
    end
    # _delete_data_object
    @testset "_delete_data_object" begin
        @test InfiniteOpt._delete_data_object(cref1) isa Nothing
        @test length(InfiniteOpt._data_dictionary(cref1)) == 1
        @test !is_valid(m, cref1)
    end
end

# Test definition methods
@testset "Definition" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_parameter(m, 0 <= pars[1:2] <= 1)
    @infinite_variable(m, inf(par))
    @infinite_variable(m, inf2(par, pars))
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, x, parameter_bounds = (par == 1))
    var = build_variable(error, inf2, Dict{Int, Float64}(1 => 0.5), check = false)
    rv = add_variable(m, var)
    data = TestData(par, 1, 1)
    meas = @measure(inf + par - x, data, name = "test")
    mindex = MeasureIndex(1)
    dinf = @deriv(inf, par)
    @testset "JuMP.build_constraint (Single)" begin
        # test bounded constraint
        bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))
        @test jump_function(build_constraint(error, inf, MOI.EqualTo(42.0),
                                             parameter_bounds = bounds)) == inf
        @test moi_set(build_constraint(error, inf, MOI.EqualTo(42.0),
                                 parameter_bounds = bounds)) == MOI.EqualTo(42.0)
        @test parameter_bounds(build_constraint(error, inf, MOI.EqualTo(42.0),
                                         parameter_bounds = bounds)) === bounds
        @test original_parameter_bounds(build_constraint(error, inf,
                       MOI.EqualTo(42.0), parameter_bounds = bounds)) !== bounds
        @test original_parameter_bounds(build_constraint(error, inf,
                        MOI.EqualTo(42.0), parameter_bounds = bounds)) == bounds
        # test scalar constraint
        con = ScalarConstraint(inf, MOI.EqualTo(42.0))
        @test jump_function(build_constraint(error, inf, MOI.EqualTo(42.0))) == inf
        @test moi_set(build_constraint(error, inf, MOI.EqualTo(42.0))) == MOI.EqualTo(42.0)
    end
    # test build_constraint (expr)
    @testset "JuMP.build_constraint (Expr)" begin
        # test bounded constraint
        bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))
        @test jump_function(build_constraint(error, inf + par, MOI.EqualTo(42.0),
                                             parameter_bounds = bounds)) == inf + par
        @test moi_set(build_constraint(error, inf + par, MOI.EqualTo(42.0),
                                 parameter_bounds = bounds)) == MOI.EqualTo(42.0)
        @test parameter_bounds(build_constraint(error, inf + par, MOI.EqualTo(42.0),
                                         parameter_bounds = bounds)) === bounds
        @test original_parameter_bounds(build_constraint(error, inf + par,
                       MOI.EqualTo(42.0), parameter_bounds = bounds)) !== bounds
        @test original_parameter_bounds(build_constraint(error, inf + par,
                        MOI.EqualTo(42.0), parameter_bounds = bounds)) == bounds
        # test scalar constraint
        con = ScalarConstraint(inf + par, MOI.EqualTo(42.0))
        @test jump_function(build_constraint(error, inf + par, MOI.EqualTo(42.0))) == inf + par
        @test moi_set(build_constraint(error, inf + par, MOI.EqualTo(42.0))) == MOI.EqualTo(42.0)
    end
    # test _check_and_update_bounds (BoundedScalarConstraint)
    @testset "_check_and_update_bounds (Bounded)" begin
        bounds = ParameterBounds()
        c = BoundedScalarConstraint(x + inf, MOI.EqualTo(0.0), bounds, copy(bounds))
        @test InfiniteOpt._check_and_update_bounds(m, c,
                                      [x, inf]).bounds[par] == IntervalSet(1, 1)
        @test original_parameter_bounds(InfiniteOpt._check_and_update_bounds(m, c,
                                                 [x, inf])) == ParameterBounds()
    end
    # test _check_and_update_bounds (ScalarConstraint)
    @testset "_check_and_update_bounds (Scalar)" begin
        c = JuMP.ScalarConstraint(inf + x, MOI.EqualTo(0.0))
        @test InfiniteOpt._check_and_update_bounds(m, c,
                            [x, inf]).bounds[par] == IntervalSet(1, 1)
        set_parameter_bounds(x, ParameterBounds(), force = true)
        m.has_hold_bounds = false
        @test InfiniteOpt._check_and_update_bounds(m, c, [x, inf]) == c
    end
    # test _update_var_constr_mapping
    @testset "_update_var_constr_mapping" begin
        # add dumby constraints
        con = ScalarConstraint(zero(AffExpr), MOI.EqualTo(0.0))
        object1 = ConstraintData(con, Int[], "c", MeasureIndex[], false)
        object2 = ConstraintData(con, Int[], "c", MeasureIndex[], false)
        idx1 = ConstraintIndex(1)
        idx2 = ConstraintIndex(2)
        cref1 = InfOptConstraintRef(m, idx1)
        cref2 = InfOptConstraintRef(m, idx2)
        @test InfiniteOpt._add_data_object(m, object1) == idx1
        @test InfiniteOpt._add_data_object(m, object2) == idx2
        # test initial use of variable
        @test isa(InfiniteOpt._update_var_constr_mapping([inf, par, meas, rv], cref1),
                  Nothing)
        @test InfiniteOpt._constraint_dependencies(inf) == [idx1]
        @test InfiniteOpt._constraint_dependencies(par) == [idx1]
        @test InfiniteOpt._constraint_dependencies(meas) == [idx1]
        @test InfiniteOpt._constraint_dependencies(rv) == [idx1]
        @test InfiniteOpt._measure_dependencies(cref1) == [mindex]
        # test secondary use of variable
        @test isa(InfiniteOpt._update_var_constr_mapping([inf, par, meas, rv], cref2),
                  Nothing)
        @test InfiniteOpt._constraint_dependencies(inf) == [idx1, idx2]
        @test InfiniteOpt._constraint_dependencies(par) == [idx1, idx2]
        @test InfiniteOpt._constraint_dependencies(meas) == [idx1, idx2]
        @test InfiniteOpt._constraint_dependencies(rv) == [idx1, idx2]
        @test InfiniteOpt._measure_dependencies(cref2) == [mindex]
        # Undo changes
        empty!(InfiniteOpt._constraint_dependencies(inf))
        empty!(InfiniteOpt._constraint_dependencies(par))
        empty!(InfiniteOpt._constraint_dependencies(meas))
        empty!(InfiniteOpt._constraint_dependencies(rv))
        empty!(InfiniteOpt._measure_dependencies(cref1))
        empty!(InfiniteOpt._measure_dependencies(cref2))
        m.constraints = MOIUC.CleverDict{ConstraintIndex, ConstraintData}()
    end
    # test add_constraint
    @testset "JuMP.add_constraint" begin
        # test reject vector constraint
        con = VectorConstraint([inf + pt], MOI.Zeros(1))
        @test_throws ErrorException add_constraint(m, con, "a")
        # test rejects bad variables
        @hold_variable(InfiniteModel(), z)
        con = ScalarConstraint(z, MOI.EqualTo(42.0))
        @test_throws VariableNotOwned{GeneralVariableRef} add_constraint(m, con, "a")
        # test bounded constraint
        bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))
        con = BoundedScalarConstraint(inf + pt, MOI.EqualTo(42.0), bounds, bounds)
        idx = ConstraintIndex(1)
        cref = InfOptConstraintRef(m, idx)
        @test add_constraint(m, con, "a") == cref
        # test infinite constraint
        con = ScalarConstraint(inf + pt, MOI.EqualTo(42.0))
        idx = ConstraintIndex(2)
        cref = InfOptConstraintRef(m, idx)
        @test add_constraint(m, con, "b") == cref
        # test finite constraint
        con = ScalarConstraint(x + pt, MOI.EqualTo(42.0))
        idx = ConstraintIndex(3)
        cref = InfOptConstraintRef(m, idx)
        @test add_constraint(m, con, "d") == cref
        @test name(cref) == "d"
        @test !InfiniteOpt._is_info_constraint(cref)
        @test !optimizer_model_ready(m)
        @test used_by_constraint(pt)
        # test with bounded hold variables
        @set_parameter_bounds(x, par == 1)
        con = ScalarConstraint(inf + pt + x, MOI.EqualTo(42.0))
        idx = ConstraintIndex(4)
        cref = InfOptConstraintRef(m, idx)
        @test add_constraint(m, con, "b") == cref
        @test parameter_bounds(cref)[par] == IntervalSet(1, 1)
        @test InfiniteOpt._core_constraint_object(cref).orig_bounds == ParameterBounds()
        set_parameter_bounds(x, ParameterBounds(), force = true)
    end
    # test macro
    @testset "JuMP.@constraint" begin
        # test scalar constraint
        idx = ConstraintIndex(5)
        cref = InfOptConstraintRef(m, idx)
        @test @constraint(m, e, x + pt -2 <= 2) == cref
        @test isa(InfiniteOpt._core_constraint_object(cref), ScalarConstraint)
        # test bounded scalar constraint
        bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))
        idx = ConstraintIndex(6)
        cref = InfOptConstraintRef(m, idx)
        @test @constraint(m, f, inf + meas - dinf <= 2, parameter_bounds = bounds) == cref
        @test InfiniteOpt._object_numbers(cref) == [1]
        @test isa(InfiniteOpt._core_constraint_object(cref), BoundedScalarConstraint)
        @test used_by_constraint(dinf)
    end
end

# Test all the methods for bound constraint macros
@testset "Bounded Constraint Macros" begin
    # test _parse_name_expression
    @testset "_parse_name_expression" begin
        @test InfiniteOpt._parse_name_expression(error, :bob) == :bob
        @test InfiniteOpt._parse_name_expression(error, :([i = 1:2])) == :([i = 1:2])
        @test InfiniteOpt._parse_name_expression(error, :(bob[i = 1:2])) == :(bob[i = 1:2])
        @test_throws ErrorException InfiniteOpt._parse_name_expression(error, :((2, 3)))
    end
    # test _extract_bounds (:call)
    @testset "_extract_bounds (:call)" begin
        # test single anonymous bound
        args = [:in, :t, :([0, 1])]
        set = :(IntervalSet(0, 1))
        tp = :(t)
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(($(dict_arg),)))
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:call)) == (nothing, bounds)
        # test with name expression
        args = [:(bob[i = 1:2]), :(t in [0, 1])]
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:call)) == (args[1], bounds)
    end
    # test @BDconstraint
    @testset "@BDconstraint" begin
        # initialize model and references
        m = InfiniteModel()
        @infinite_parameter(m, 0 <= par <= 10)
        @infinite_parameter(m, 0 <= pars[1:2] <= 10)
        @infinite_variable(m, inf(par))
        @point_variable(m, inf(0.5), pt)
        @hold_variable(m, x)
        bounds1 = ParameterBounds(Dict(par => IntervalSet(0, 1)))
        m2 = Model()
        # test errors
        @test_macro_throws ErrorException @BDconstraint(m, par = 0, inf == 0,
                                                     parameter_bounds = bounds1)
        @test_macro_throws ErrorException @BDconstraint(m, inf == 0)
        @test_macro_throws ErrorException @BDconstraint(m, con, inf == 0)
        @test_macro_throws ErrorException @BDconstraint(m, con[1:2], inf == 0)
        @test_macro_throws ErrorException @BDconstraint(m, [1:2], inf == 0)
        @test_macro_throws ErrorException @BDconstraint(m, con, inf == 0)
        @test_macro_throws ErrorException @BDconstraint(m, a.b, inf == 0)
        @test_macro_throws ErrorException @BDconstraint(m2, par == 0, inf == 0)
        @test_macro_throws ErrorException @BDconstraint(m2, con(par == 0), inf == 0)
        # test anonymous constraint with set
        idx = ConstraintIndex(1)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, par in [0, 1], inf + x == 2) == cref
        @test parameter_bounds(cref) == bounds1
        # test anonymous constraint with comparison
        idx = ConstraintIndex(2)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, 0 <= par <= 1, inf + x == 2) == cref
        @test parameter_bounds(cref) == bounds1
        # test anonymous constraint with equality
        idx = ConstraintIndex(3)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, par == 0, inf + x == 2) == cref
        @test parameter_bounds(cref) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test anonymous multiple constraint
        idxs = [ConstraintIndex(4), ConstraintIndex(5)]
        crefs = [InfOptConstraintRef(m, idx) for idx in idxs]
        @test @BDconstraint(m, [1:2](par == 0), inf + x == 2) == crefs
        @test parameter_bounds(crefs[1]) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test parameter_bounds(crefs[2]) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test anonymous multiple bounds
        idx = ConstraintIndex(6)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, (par == 0, pars[1] in [0, 1]), inf + x == 2) == cref
        @test parameter_bounds(cref) == ParameterBounds((par => IntervalSet(0, 0),
                                                  pars[1] => IntervalSet(0, 1)))
        # test anonymous multiple bounds with vector input
        idx = ConstraintIndex(7)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, (par == 0, pars in [0, 1]), inf + x == 2) == cref
        @test parameter_bounds(cref) == ParameterBounds((par => IntervalSet(0, 0),
                    pars[1] => IntervalSet(0, 1), pars[2] => IntervalSet(0, 1)))
        # test named constraint
        idx = ConstraintIndex(8)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, c1(par in [0, 1]), inf + x == 2) == cref
        @test parameter_bounds(cref) == bounds1
        # test named constraint with other vector type
        idx = ConstraintIndex(9)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, c2(par in [0, 1], 0 <= pars <= 1), inf + x == 2) == cref
        @test parameter_bounds(cref) == ParameterBounds((par => IntervalSet(0, 1),
                    pars[1] => IntervalSet(0, 1), pars[2] => IntervalSet(0, 1)))
        # test named constraints
        idxs = [ConstraintIndex(10), ConstraintIndex(11)]
        crefs = [InfOptConstraintRef(m, idx) for idx in idxs]
        @test @BDconstraint(m, c3[1:2](par == 0), inf + x == 2) == crefs
        @test parameter_bounds(crefs[1]) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test parameter_bounds(crefs[2]) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test different container type
        idxs = [ConstraintIndex(12), ConstraintIndex(13)]
        crefs = [InfOptConstraintRef(m, idx) for idx in idxs]
        expected = convert(JuMPC.SparseAxisArray, crefs)
        @test @BDconstraint(m, c4[1:2](par == 0), inf + x == 2,
                            container = SparseAxisArray) == expected
        @test parameter_bounds(crefs[1]) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test parameter_bounds(crefs[2]) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test with SparseAxisArray input
        idx = ConstraintIndex(14)
        cref = InfOptConstraintRef(m, idx)
        pars = convert(JuMPC.SparseAxisArray, pars)
        @test @BDconstraint(m, c5(pars in [0, 1]), inf + x == 2) == cref
        @test parameter_bounds(cref) == ParameterBounds((pars[1] => IntervalSet(0, 1),
                                                  pars[2] => IntervalSet(0, 1)))
        # test with hold variable bounds
        @set_parameter_bounds(x, pars[1] in [0, 2])
        idx = ConstraintIndex(15)
        cref = InfOptConstraintRef(m, idx)
        @test @BDconstraint(m, c10(par in [0, 1]), inf + x == 2) == cref
        @test parameter_bounds(cref) == ParameterBounds(Dict(par => IntervalSet(0, 1),
                                                  pars[1] => IntervalSet(0, 2)))
        @test InfiniteOpt._core_constraint_object(cref).orig_bounds == ParameterBounds(Dict(par => IntervalSet(0, 1)))
    end
end

# Test parameter reference methods
@testset "Parameter References" begin
    m = InfiniteModel()
    @independent_parameter(m, t in [0, 1])
    @independent_parameter(m, y in [0, 1])
    @dependent_parameters(m, x[1:3] in [0, 1])
    @hold_variable(m, z)
    @constraint(m, c1, z >= 0)
    @constraint(m, c2, z + t + x[1] >= 0)
    # test parameter_refs
    @testset "parameter_refs" begin
        @test parameter_refs(c1) == ()
        @test parameter_refs(c2) == (t, x)
    end
end

# Test parameter bound methods
@testset "Parameter Bounds" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 10)
    @infinite_parameter(m, 0 <= pars[1:2] <= 10)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, x >= 0, Int)
    @BDconstraint(m, c1(par in [0, 1]), inf + x == 0)
    @constraint(m, c2, x * pt + x == 2)
    # test _update_constr_param_bounds
    @testset "_update_constr_param_bounds" begin
        # test empty bounds
        bounds = ParameterBounds()
        @test isa(InfiniteOpt._update_constr_param_bounds(c1, bounds, bounds),
                  Nothing)
        @test isa(constraint_object(c1), ScalarConstraint)
        # test normalized_coefficient
        bounds = ParameterBounds(Dict(par => IntervalSet(0, 5)))
        @test isa(InfiniteOpt._update_constr_param_bounds(c1, bounds, bounds),
                  Nothing)
        @test parameter_bounds(c1) == bounds
    end
    # test set_parameter_bounds
    @testset "set_parameter_bounds" begin
        # test force error
        bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))
        @test_throws ErrorException set_parameter_bounds(c1, bounds)
        # test normal
        @test isa(set_parameter_bounds(c2, bounds), Nothing)
        @test parameter_bounds(c2) == bounds
        @test !optimizer_model_ready(m)
        # test test error with bounds
        bounds = ParameterBounds(Dict(par => IntervalSet(-1, 1)))
        @test_throws ErrorException set_parameter_bounds(c1, bounds, force = true)
    end
    # test add_parameter_bound
    @testset "add_parameter_bounds" begin
        # test already has bounds
        bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))
        @test isa(add_parameter_bounds(c1, copy(bounds)), Nothing)
        @test parameter_bounds(c1) == bounds
        # test doesn't have bounds
        @constraint(m, c3, inf + pt == 0)
        @test isa(add_parameter_bounds(c3, copy(bounds)), Nothing)
        @test parameter_bounds(c3) == bounds
    end
    # test delete_parameter_bounds
    @testset "delete_parameter_bounds" begin
        # test partial deletion
        bounds = ParameterBounds(Dict(pars[1] => IntervalSet(0, 1)))
        @test isa(add_parameter_bounds(c1, bounds), Nothing)
        @add_parameter_bounds(x, par == 0)
        @test isa(delete_parameter_bounds(c1), Nothing)
        expected = ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test parameter_bounds(c1) == expected
        # test again without hold bounds
        delete_parameter_bounds(x)
        @test isa(constraint_object(c1), ScalarConstraint)
        # test already gone
        @test isa(delete_parameter_bounds(c1), Nothing)
        @test isa(constraint_object(c1), ScalarConstraint)
    end
    # test @set_parameter_bounds
    @testset "@set_parameter_bounds" begin
        # Note errors were already checked with hold variables
        # test force error
        @test_macro_throws ErrorException @set_parameter_bounds(c1, par == 1)
        # test with single
        @set_parameter_bounds(x, pars == 0)
        @test isa(@set_parameter_bounds(c1, par == 1, force = true), Nothing)
        @test parameter_bounds(c1) == ParameterBounds((par => IntervalSet(1, 1),
                                                  pars => IntervalSet(0, 0)))
        set_parameter_bounds(x, ParameterBounds(), force = true)
        # test multiple
        @test isa(@set_parameter_bounds(c1, pars == 0, force = true), Nothing)
        @test parameter_bounds(c1)[pars[2]] == IntervalSet(0, 0)
    end
    # test @add_parameter_bounds
    @testset "@add_parameter_bounds" begin
        # Note errors were already checked with hold variables
        # test bounds error
        @test_macro_throws ErrorException @add_parameter_bounds(c3, par in [-1, 1])
        # test with multiple
        @test isa(@add_parameter_bounds(c2, (pars == 0, par in [0, 1])), Nothing)
        @test parameter_bounds(c2)[pars[2]] == IntervalSet(0, 0)
        @test parameter_bounds(c2)[par] == IntervalSet(0, 1)
    end
end

# Test coefficient queries and modifications
@testset "Coefficient Methods" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, x >= 0, Int)
    @BDconstraint(m, c1(par in [0, 1]), inf + x == 0)
    @constraint(m, c2, x * pt + x == 2)
    # test set_normalized_coefficient
    @testset "JuMP.set_normalized_coefficient" begin
        @test isa(set_normalized_coefficient(c1, x, 2), Nothing)
        @test isa(constraint_object(c1), BoundedScalarConstraint)
        @test jump_function(constraint_object(c1)) == inf + 2x
        @test isa(set_normalized_coefficient(c2, inf, 2), Nothing)
        @test isa(constraint_object(c2), ScalarConstraint)
        @test jump_function(constraint_object(c2)) == x * pt + x + 2inf
    end
    # test normalized_coefficient
    @testset "JuMP.normalized_coefficient" begin
        @test normalized_coefficient(c1, x) == 2.
        @test normalized_coefficient(c1, pt) == 0.
        @test normalized_coefficient(c2, inf) == 2.
        @test normalized_coefficient(LowerBoundRef(x), x) == 1.
        @test normalized_coefficient(LowerBoundRef(x), pt) == 0.
    end
    # test _set_set_value
    @testset "_set_set_value" begin
        @test InfiniteOpt._set_set_value(MOI.EqualTo(1.), 2) == MOI.EqualTo(2.)
        @test InfiniteOpt._set_set_value(MOI.GreaterThan(1.), 2) == MOI.GreaterThan(2.)
        @test InfiniteOpt._set_set_value(MOI.LessThan(1.), 2) == MOI.LessThan(2.)
    end
    # test set_normalized_rhs
    @testset "JuMP.set_normalized_rhs" begin
        @test isa(set_normalized_rhs(c1, 42), Nothing)
        @test isa(constraint_object(c1), BoundedScalarConstraint)
        @test constraint_object(c1).set == MOI.EqualTo(42.0)
        @test isa(set_normalized_rhs(c2, 42), Nothing)
        @test isa(constraint_object(c2), ScalarConstraint)
        @test constraint_object(c2).set == MOI.EqualTo(42.0)
        @test isa(set_normalized_rhs(LowerBoundRef(x), 42), Nothing)
        @test isa(constraint_object(LowerBoundRef(x)), ScalarConstraint)
        @test constraint_object(LowerBoundRef(x)).set == MOI.GreaterThan(42.0)
    end
    # test normalized_rhs
    @testset "JuMP.normalized_rhs" begin
        @test normalized_rhs(c1) == 42.
        @test normalized_rhs(c2) == 42.
        @test normalized_rhs(LowerBoundRef(x)) == 42.
    end
    # test add_to_function_constant
    @testset "JuMP.add_to_function_constant" begin
        @test isa(add_to_function_constant(c1, 4), Nothing)
        @test isa(constraint_object(c1), BoundedScalarConstraint)
        @test constraint_object(c1).set == MOI.EqualTo(38.0)
        @test isa(add_to_function_constant(c2, 4), Nothing)
        @test isa(constraint_object(c2), ScalarConstraint)
        @test constraint_object(c2).set == MOI.EqualTo(38.0)
        @test isa(add_to_function_constant(LowerBoundRef(x), 4), Nothing)
        @test isa(constraint_object(LowerBoundRef(x)), ScalarConstraint)
        @test constraint_object(LowerBoundRef(x)).set == MOI.GreaterThan(38.0)
    end
end

# Test num_constraints extensions
@testset "JuMP.num_constraints Extensions" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par) >= 0)
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, 0 <= x <= 1, Int)
    # test search function type and set type
    @testset "Function and Set" begin
        @test num_constraints(m, GeneralVariableRef, MOI.LessThan) == 1
        @test num_constraints(m, GeneralVariableRef, MOI.GreaterThan) == 3
    end
    # test search function type
    @testset "Function" begin
        @test num_constraints(m, GeneralVariableRef) == 5
    end
    # test search set type
    @testset "Set" begin
        @test num_constraints(m, MOI.LessThan) == 1
        @test num_constraints(m, MOI.GreaterThan) == 3
    end
    # test search total
    @testset "Total" begin
        @test num_constraints(m) == 5
    end
end

# Test all_constraints extensions
@testset "JuMP.all_constraints Extensions" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par) >= 0)
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, 0 <= x <= 1, Int)
    # test search function type and set type
    @testset "Function and Set" begin
        list = [InfOptConstraintRef(m, ConstraintIndex(4))]
        @test all_constraints(m, GeneralVariableRef, MOI.LessThan) == list
        list = [InfOptConstraintRef(m, ConstraintIndex(1)),
                InfOptConstraintRef(m, ConstraintIndex(2)),
                InfOptConstraintRef(m, ConstraintIndex(3))]
        @test all_constraints(m, GeneralVariableRef, MOI.GreaterThan) == list
    end
    # test search function type
    @testset "Function" begin
        list = [InfOptConstraintRef(m, ConstraintIndex(1)),
                InfOptConstraintRef(m, ConstraintIndex(2)),
                InfOptConstraintRef(m, ConstraintIndex(3)),
                InfOptConstraintRef(m, ConstraintIndex(4)),
                InfOptConstraintRef(m, ConstraintIndex(5))]
        @test all_constraints(m, GeneralVariableRef) == list
    end
    # test search set type
    @testset "Set" begin
        list = [InfOptConstraintRef(m, ConstraintIndex(4))]
        @test all_constraints(m, MOI.LessThan) == list
        list = [InfOptConstraintRef(m, ConstraintIndex(1)),
                InfOptConstraintRef(m, ConstraintIndex(2)),
                InfOptConstraintRef(m, ConstraintIndex(3))]
        @test all_constraints(m, MOI.GreaterThan) == list
    end
    # test search total
    @testset "Total" begin
        list = [InfOptConstraintRef(m, ConstraintIndex(1)),
                InfOptConstraintRef(m, ConstraintIndex(2)),
                InfOptConstraintRef(m, ConstraintIndex(3)),
                InfOptConstraintRef(m, ConstraintIndex(4)),
                InfOptConstraintRef(m, ConstraintIndex(5))]
        @test all_constraints(m) == list
    end
end

# Test list_of_constraint_types
@testset "JuMP.list_of_constraint_types" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par) >= 0)
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, 0 <= x <= 1, Int)
    expected = [(GeneralVariableRef, MOI.GreaterThan{Float64}),
                (GeneralVariableRef, MOI.LessThan{Float64}),
                (GeneralVariableRef, MOI.Integer)]
    @test isempty(setdiff(list_of_constraint_types(m), expected))
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

prime_movers = ["recip_engine", "micro_turbine", "combustion_turbine", "fuel_cell"]

"""
    CHP

struct with outer constructor:
```julia
prime_mover::String = ""
# following must be provided by user if not providing prime_mover
installed_cost_per_kw::Union{Float64, AbstractVector{Float64}} = NaN
tech_sizes_for_cost_curve::Union{Float64, AbstractVector{Float64}} = NaN
om_cost_per_kwh::Float64 = NaN
elec_effic_half_load = NaN
elec_effic_full_load::Float64 = NaN
min_turn_down_pct::Float64 = NaN
thermal_effic_full_load::Float64 = NaN
thermal_effic_half_load::Float64 = NaN
min_allowable_kw::Float64 = NaN
max_kw::Float64 = NaN
cooling_thermal_factor::Float64 = NaN  # only needed with cooling load
unavailability_periods::AbstractVector{Dict} = Dict[]

size_class::Int = 1
min_kw::Float64 = 0.0
fuel_cost_per_mmbtu::Union{<:Real, AbstractVector{<:Real}} = 0.0,
om_cost_per_kw::Float64 = 0.0
om_cost_per_hr_per_kw_rated::Float64 = 0.0
supplementary_firing_capital_cost_per_kw::Float64 = 150.0
supplementary_firing_max_steam_ratio::Float64 = 1.0
supplementary_firing_efficiency::Float64 = 0.92
standby_rate_us_dollars_per_kw_per_month = 0.0
reduces_demand_charges = true
use_default_derate::Bool = true
max_derate_factor::Float64 = 1.0
derate_start_temp_degF::Float64 = 0.0
derate_slope_pct_per_degF::Float64 = 0.0
can_supply_steam_turbine::Bool=false

macrs_option_years::Int = 5
macrs_bonus_pct::Float64 = 1.0
macrs_itc_reduction::Float64 = 0.5
federal_itc_pct::Float64 = 0.1
federal_rebate_per_kw::Float64 = 0.0
state_ibi_pct::Float64 = 0.0
state_ibi_max::Float64 = 1.0e10
state_rebate_per_kw::Float64 = 0.0
state_rebate_max::Float64 = 1.0e10
utility_ibi_pct::Float64 = 0.0
utility_ibi_max::Float64 = 1.0e10
utility_rebate_per_kw::Float64 = 0.0
utility_rebate_max::Float64 = 1.0e10
production_incentive_per_kwh::Float64 = 0.0
production_incentive_max_benefit::Float64 = 1.0e9
production_incentive_years::Int = 0
production_incentive_max_kw::Float64 = 1.0e9
can_net_meter::Bool = false
can_wholesale::Bool = false
can_export_beyond_nem_limit::Bool = false
can_curtail::Bool = false
# emissions_factor_lb_CO2_per_mmbtu::Float64
```
"""
Base.@kwdef mutable struct CHP <: AbstractCHP
    prime_mover::String = ""
    # following must be provided by user if not providing prime_mover
    installed_cost_per_kw::Union{Float64, AbstractVector{Float64}} = Float64[]
    tech_sizes_for_cost_curve::AbstractVector{Float64} = Float64[]
    om_cost_per_kwh::Float64 = NaN
    elec_effic_half_load = NaN
    elec_effic_full_load::Float64 = NaN
    min_turn_down_pct::Float64 = NaN
    thermal_effic_full_load::Float64 = NaN
    thermal_effic_half_load::Float64 = NaN
    min_allowable_kw::Float64 = NaN
    max_kw::Float64 = NaN
    cooling_thermal_factor::Float64 = NaN  # only needed with cooling load
    unavailability_periods::AbstractVector{Dict} = Dict[]

    size_class::Int = 1
    min_kw::Float64 = 0.0
    fuel_cost_per_mmbtu::Union{<:Real, AbstractVector{<:Real}} = 0.0
    om_cost_per_kw::Float64 = 0.0
    om_cost_per_hr_per_kw_rated::Float64 = 0.0
    supplementary_firing_capital_cost_per_kw::Float64 = 150.0
    supplementary_firing_max_steam_ratio::Float64 = 1.0
    supplementary_firing_efficiency::Float64 = 0.92
    standby_rate_us_dollars_per_kw_per_month::Float64 = 0.0
    reduces_demand_charges::Bool = true
    use_default_derate::Bool = true
    max_derate_factor::Float64 = 1.0
    derate_start_temp_degF::Float64 = 0.0
    derate_slope_pct_per_degF::Float64 = 0.0
    can_supply_steam_turbine::Bool=false

    macrs_option_years::Int = 5
    macrs_bonus_pct::Float64 = 1.0
    macrs_itc_reduction::Float64 = 0.5
    federal_itc_pct::Float64 = 0.1
    federal_rebate_per_kw::Float64 = 0.0
    state_ibi_pct::Float64 = 0.0
    state_ibi_max::Float64 = 1.0e10
    state_rebate_per_kw::Float64 = 0.0
    state_rebate_max::Float64 = 1.0e10
    utility_ibi_pct::Float64 = 0.0
    utility_ibi_max::Float64 = 1.0e10
    utility_rebate_per_kw::Float64 = 0.0
    utility_rebate_max::Float64 = 1.0e10
    production_incentive_per_kwh::Float64 = 0.0
    production_incentive_max_benefit::Float64 = 1.0e9
    production_incentive_years::Int = 0
    production_incentive_max_kw::Float64 = 1.0e9
    can_net_meter::Bool = false
    can_wholesale::Bool = false
    can_export_beyond_nem_limit::Bool = false
    can_curtail::Bool = false
    # emissions_factor_lb_CO2_per_mmbtu::Float64
end

function CHP(d::Dict)
    chp = CHP(; dictkeys_tosymbols(d)...)

    # Must provide prime_mover or all of custom_chp_inputs
    custom_chp_inputs = Dict{Symbol, Any}(
        :installed_cost_per_kw => chp.installed_cost_per_kw,
        :tech_sizes_for_cost_curve => chp.tech_sizes_for_cost_curve,
        :om_cost_per_kwh => chp.om_cost_per_kwh,
        :elec_effic_full_load => chp.elec_effic_full_load,
        :elec_effic_half_load => chp.elec_effic_half_load,
        :min_turn_down_pct => chp.min_turn_down_pct,
        :thermal_effic_full_load => chp.thermal_effic_full_load,
        :thermal_effic_half_load => chp.thermal_effic_half_load,
        :min_allowable_kw => chp.min_allowable_kw,
        :max_kw => chp.max_kw,
        :cooling_thermal_factor => chp.cooling_thermal_factor
    )

    # Installed cost input validation
    update_installed_cost_params = false
    pass_all_params_error = false
    if !isempty(chp.installed_cost_per_kw) && typeof(chp.installed_cost_per_kw) == Float64
        if !isempty(chp.tech_sizes_for_cost_curve)
            chp.tech_sizes_for_cost_curve = []
            @warn "Ignoring `chp.tech_sizes_for_cost_curve` input because `chp.installed_cost_per_kw` is a scalar"
        end
    elseif !isempty(chp.installed_cost_per_kw) && isempty(chp.tech_sizes_for_cost_curve)
        @error "To model CHP cost curve, you must provide `chp.tech_sizes_for_cost` vector of equal length to `chp.installed_cost_per_kw`"
    elseif isempty(chp.tech_sizes_for_cost_curve)
        update_installed_cost_params = true
    elseif isempty(chp.prime_mover)
        pass_all_params_error = true
    end

    if isempty(chp.prime_mover)
        if !pass_all_params_error
            if any(isnan(v) for v in values(custom_chp_inputs)) || isempty(chp.unavailability_periods)
                pass_all_params_error = true
            end
        end
        if pass_all_params_error
            @error "To model CHP you must provide at least `prime_mover` from $(prime_movers) or all of $([string(k) for k in keys(custom_chp_inputs)]) and unavailability_periods."
        end
    elseif !(isempty(chp.prime_mover))
        @assert chp.prime_mover in prime_movers
        if chp.prime_mover == "combustion_turbine"
            boiler_type = "steam"
        else
            boiler_type = "hot_water"
        end
        # set all missing default values in custom_chp_inputs
        defaults = get_prime_mover_defaults(chp.prime_mover, boiler_type, chp.size_class)
        for (k, v) in custom_chp_inputs
            if k in [:installed_cost_per_kw, :tech_sizes_for_cost_curve]
                if update_installed_cost_params
                    setproperty!(chp, k, defaults[string(k)])
                end
            elseif isnan(v)
                setproperty!(chp, k, defaults[string(k)])
            end
        end
        if isempty(chp.unavailability_periods)
            chp.unavailability_periods = defaults["unavailability_periods"]
        end
    end

    return chp
end

"""
    get_prime_mover_defaults(prime_mover::String, boiler_type::String, size_class::Int)

return a Dict{String, Union{Float64, AbstractVector{Float64}}} by selecting the appropriate values from
data/chp/chp_default_data.json, which contains values based on prime_mover, boiler_type, and size_class for the
custom_chp_inputs, i.e.
- "installed_cost_per_kw"
- "tech_sizes_for_cost_curve"
- "om_cost_per_kwh"
- "elec_effic_full_load"
- "min_turn_down_pct",
- "thermal_effic_full_load"
- "thermal_effic_half_load"
- "unavailability_periods"
"""
function get_prime_mover_defaults(prime_mover::String, boiler_type::String, size_class::Int)
    pmds = JSON.parsefile(joinpath(dirname(@__FILE__), "..", "..", "data", "chp", "chp_defaults.json"))
    prime_mover_defaults = Dict{String, Any}()

    for key in keys(pmds[prime_mover])
        if key in ["thermal_effic_full_load", "thermal_effic_half_load"]
            prime_mover_defaults[key] = pmds[prime_mover][key][boiler_type][size_class]
        elseif key == "unavailability_periods"
            prime_mover_defaults[key] = convert(Vector{Dict}, pmds[prime_mover][key])
        else
            prime_mover_defaults[key] = pmds[prime_mover][key][size_class]
        end
    end
    pmds = nothing

    for (k,v) in prime_mover_defaults
        if typeof(v) <: AbstractVector{Any} && k != "unavailability_periods"
            prime_mover_defaults[k] = convert(Vector{Float64}, v)  # JSON.parsefile makes things Vector{Any}
        end
    end
    return prime_mover_defaults
end

ext: jl

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
# Adapted from https://github.com/arviz-devs/arviz/blob/master/doc/logo/generate_logo.py
# Originally licensed Under Apache-2.0 to ArviZ developers.
using Distributions, PyCall, PyPlot

# Config options
prefix = "logo"
extensions = ["png", "pdf"]
dist_colors = ["#72D0F5", "#72D0F5", "#4063D8"] # right to left
dot_colors = ["#cb3c33", "#9558b2", "#389826"] # counter-clockwise
dot_coords = [0.3065 0.764; 0.3415 0.764; 0.32425 1.012]
dot_radius = 120

# Setup
bbox = matplotlib.transforms.Bbox([0.75 0.5; 5.4 2.2])
save_kwargs = (dpi = 300, bbox_inches = bbox, transparent = true)
cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", dist_colors)

axis("off")
ylim(0, 5.5)
xlim(0, 0.9)

# Draw distribution
x = range(0, 1; length = 200)
pdfx = pdf.(Beta(2, 5), x)

path = matplotlib.path.Path([x pdfx])
patch = matplotlib.patches.PathPatch(path; facecolor = "none", alpha = 0)
gca().add_patch(patch)

imshow(
    [1 0 0; 1 1 0];
    cmap = cmap,
    interpolation = "bicubic",
    origin = "lower",
    extent = [0, 1, 0.0, 5],
    aspect = "auto",
    clip_path = patch,
    clip_on = true,
    zorder = 0,
)

# Save text-free logo
savefig("$(prefix)_notext.png"; save_kwargs...)

# Add text
logotext = text(
    x = 0.04,
    y = -0.01,
    s = "ArviZ",
    clip_on = true,
    fontdict = Dict("name" => "ubuntu mono", "fontsize" => 62),
    color = "w",
    zorder = 1,
)

# Add Julia dots
for i = 1:3
    scatter(dot_coords[i, :]...; s = dot_radius, zorder = 2, color = dot_colors[i])
end

# Save logo
for ext in extensions
    savefig("$(prefix).$(ext)"; save_kwargs...)
end

logotext.set_color("black")
for ext in extensions
    savefig("$(prefix)-dark.$(ext)"; save_kwargs...)
end


text:
############################################################
# Linked list with last (without nothing)
############################################################

mutable struct Orbit3{T, S} <: AbstractOrbit{T,S}
    pt::T # starting point
    last::Tuple{T, S} # idx, content
    elts::Dict{T, Tuple{S, T}} # idx -> (content, nextptr)
end

function Orbit3(pt::T, val::S) where {T,S}
    return Orbit3(pt, (pt, val), Dict{T, Tuple{S, T}}())
end

Orbit3(pt::T) where {T} = Orbit3(pt, nothing)

@inline Base.in(x, orb::Orbit3) = haskey(orb.elts, x) || x == orb.last[1]

@inline function Base.push!(orb::Orbit3{T,S}, tup::Tuple{T,S}) where {T,S}
    new, val = tup
    orb.elts[first(orb.last)] = (last(orb.last), new)
    orb.last = (new, val)
    return orb
end

@inline Base.push!(orb::Orbit3{T,Nothing}, n::T) where T = push!(orb, (n,nothing))

# o is a point in orbit, return (value, nextptr)
@inline Base.getindex(orb::Orbit3, o) = (islast(orb, o) ? last(orb.last) : first(orb.elts[o]))

@inline islast(orb::Orbit3, s) = s == first(orb.last)
@inline next(orb::Orbit3, s) = (last(orb.elts[s]), last(orb.elts[s]))

@inline Base.iterate(orb::Orbit3) = orb.pt, orb.pt
@inline Base.iterate(orb::Orbit3, s) = (islast(orb, s) ? nothing : next(orb, s))
@inline Base.length(orb::Orbit3) = length(orb.elts) + 1

@inline function Base.:(==)(o1::Orbit3, o2::Orbit3)
    o1.pt == o2.pt && o1.last == o2.last || return false
    return o1.elts == o2.elts
end

@inline Base.first(orb::Orbit3) = orb.pt


text:
using Richardson, Test, LinearAlgebra

@testset "basic tests" begin
    X = Float64[]
    val, err = extrapolate(1.0, rtol=1e-10, contract=0.1) do x
        push!(X, x)
        sin(x)/x
    end
    @test val ≈ 1 rtol=1e-10
    @test err < 1e-10
    @test all(x -> x > 0, X)
    @test length(X) == 6

    X2 = Float64[]
    val, err = extrapolate(1.0, rtol=1e-10, contract=0.1, maxeval=3) do x
        push!(X2, x)
        sin(x)/x
    end
    @test X2 == X[1:3]
    @test val ≈ 1 rtol=1e-2
    @test err < 1e-2

    empty!(X2)
    val, err = extrapolate(1.0, rtol=1e-13, contract=0.1, power=2) do x
        push!(X2, x)
        sin(x)/x
    end
    @test X2 == X[1:5]
    @test val ≈ 1 rtol=1e-13
    @test err < 1e-13

    empty!(X)
    val, err = extrapolate(-1.0, rtol=1e-10, contract=0.1) do x
        push!(X, x)
        sin(x)/x
    end
    @test val ≈ 1 rtol=1e-10
    @test err < 1e-10
    @test all(x -> x < 0, X)
    @test length(X) == 6

    empty!(X)
    val, err = extrapolate(1.0, x0=Inf) do x
        push!(X, x)
        (x^2 + 3x - 2) / (x^2 + 5)
    end
    @test val ≈ 1 rtol=1e-12
    @test err < 1e-10
    @test length(X) == 7

    # make sure this terminates rather than looping on NaN values
    val, err = extrapolate(x -> log(x)/sqrt(x-1), 0.2, x0=1.0)
    @test abs(val) < 1e-6 && abs(err) < 1e-6

    # rapid convergence via Puiseux series:
    val, err = extrapolate(x -> log(x)/sqrt(x-1), 0.2, x0=1.0, power=0.5)
    @test abs(val) < 1e-12 && abs(err) < 1e-12

    # vector-valued function support
    val, err = extrapolate(x -> [sin(x)/x, cos(x)], 0.1, rtol=1e-10)
    @test val ≈ [1,1] rtol=1e-10
    @test err < norm([1,1])*1e-10

    # Float32 specified via starting h
    val, err = extrapolate(x -> sin(x)/x, 1.0f0, rtol=1e-5)
    @test val ≈ 1.0 rtol=1e-5
    @test err < 1e-7
    @test val isa Float32
    @test err isa Float32

    # extrapolation in a vector space
    val, err = extrapolate([0.1,-0.2],x0=[0.0,1.0]) do x
        sin(x[1])/x[1]*log(x[2])/(x[2]-1)
    end
    @test val ≈ 1 rtol=1e-12
    @test err < 1e-8

    # extrapolation for complex arguments
    val, err = extrapolate(sin, 1+2im, x0=4+5im)
    @test err < abs(val)*1e-8
    @test val ≈ sin(4+5im) rtol=1e-12

    # extrapolation for real arguments with complex contraction (spirals inwards)
    val, err = extrapolate(sin, 1, x0=4, contract = 0.5 + 0.5im)
    @test err < abs(val)*1e-8
    @test val ≈ sin(4) rtol=1e-10

    # compare to some known results
    val, err = extrapolate(x -> sin(x)/x, 1, contract=0.5, maxeval=4)
    @test 1 - val ≈ 0.0001228230317512624
    @test err ≈ 0.001886562265816072
    val, err = extrapolate(x -> sin(x)/x, 1, contract=0.5, maxeval=4, power=2)
    @test 1 - val ≈ 6.647181693963944e-10
    @test err ≈ 3.0435620007596143e-6

    # increase breaktol to fix convergence for functions with poles:
    @test -10000 ≈ extrapolate(1000.0, breaktol=Inf) do h
        x = 0.01
        (1/(x+h) - 1/x) / h
    end[1]
end

@testset "iterator API" begin
    # collection of irregular contraction ratios:
    val, err = extrapolate((sin(h)/h,h) for h in [1, 0.8, 0.3, 0.2, 0.1, 0.03, 0.01])
    @test err < 1e-7
    @test val - 1 ≈ 9.49995637711254e-12

    # non-decreasing |h| is an error:
    @test_throws ArgumentError extrapolate((sin(h)/h,h) for h in [1, 0.8, 1.3])
    @test_throws ArgumentError extrapolate((sin(h)/h,h) for h in [1, 0.8, 0.8])
end

@testset "in-place API" begin
    # collection of irregular contraction ratios:
    val, err = extrapolate!([(sin(h)/h,h) for h in [1, 0.8, 0.3, 0.2, 0.1, 0.03, 0.01]])
    @test err < 1e-7
    @test val - 1 ≈ 9.49995637711254e-12

    # non-decreasing |h| is an error:
    @test_throws ArgumentError extrapolate!([(sin(h)/h,h) for h in [1, 0.8, 1.3]])
    @test_throws ArgumentError extrapolate!([(sin(h)/h,h) for h in [1, 0.8, 0.8]])
end


text:
"""
Description of individual model.
"""
module individual
importall common
using fitness, core, utility
export sort, sort!, isbetter, best!, best

# Load contents.
include("individual/stage.jl")

"""
Sorts a list of individuals in descending order of fitness, from best to worst.
"""
sort{F}(s::FitnessScheme, l::Vector{Tuple{Int, F}}) =
  sort(l, lt = (x, y) -> isbetter(s, x, y))

"""
Sorts a list of individuals in descending order of fitness, from best to worst,
in place.
"""
sort!{F}(s::FitnessScheme, l::Vector{Tuple{Int, F}}) =
  sort!(l, lt = (x, y) -> isbetter(s, x, y))

"""
Compares an individual X to an individual Y to determine whether X is an
improvement upon Y.
"""
isbetter{F}(s::FitnessScheme, x::Tuple{Int, F}, y::Tuple{Int, F}) =
  compare(s, x, y) == -1
compare{F}(s::FitnessScheme, x::Tuple{Int, F},  y::Tuple{Int, F}) =
  compare(s, x[2], y[2])

"""
Returns the "best" individual from a provided list, according to a given
fitness scheme.

Fails if the if list of individuals is empty.
"""
function best{F}(s::FitnessScheme, inds::Vector{Tuple{Int, F}})
  bst = inds[1]
  for i in inds
    if isbetter(s, i, bst)
      bst = i
    end
  end
  return bst
end

"""
Performs any necessary post-processing on the individuals using this fitness
scheme.
"""
scale!{F}(::FitnessScheme, ::Vector{Tuple{Int, F}}) = return
end


text:
module OptMapsNLopt

using DocStringExtensions
using Random
using Reexport
using Base.Iterators
using ColorSchemes
using ColorTypes
using Images
using Luxor
using CSV
using DataFrames
@reexport using MultivariateStats
@reexport using NLopt
@reexport using JuMP

export OptMapNLopt, create_map!, draw_solution!, save_all, save_map,
       save_base_img

"""
$(TYPEDEF)

Defines the optimizer structure for NLopt solved problems

$(TYPEDFIELDS)

"""
mutable struct OptMapNLopt
    variable_names::Array{String,1}
    variable_domains::Array{Tuple{R,R},1} where R <: Number
    constraint_names::Array{String,1}
    objective_f::Function
    objective_f_grad
    objective_f_args
    constraints_f
    constraints_f_grads
    constraints_f_args
    sample_points_budget::Int
    verbose::Bool
    min_2d::AbstractArray{Float64,1}
    max_2d::AbstractArray{Float64,1}
    X_nd::AbstractArray{Float64,2}
    C::AbstractArray{Float64,2}
    Z::AbstractArray{Float64,2}
    M::Array{RGB{Normed{UInt8,8}},2}
    S::Array{RGB{Normed{UInt8,8}},2}
    pca #::MultivariateStats.PCA{Float64}
    MAP_SIDE::Int64
end

"""
$(SIGNATURES)

Constructor for the OptMapNLopt structure

Arguments:
=========

variable_names::Array{String, 1}:               array of strings naming the
                                                variables
variable_domains::Array{Tuple{R, R}, 1}         [[a, b]] for all variables
where R <: Numeber:
objective_f::Function:                          objective function
constraint_names:                               [(x) -> true]
sample_points_budget:                           number points to sample for each
                                                variable
verbose:                                        generate more or less verbose
                                                information

Returns:
=======

O::OptMapNLopt:                                 struct object for the
                                                optimization map

"""
function OptMapNLopt(variable_names::Array{String,1},
                     variable_domains::Array{Tuple{R,R},1} where R <: Number,
                     objective_f::Function;
                     objective_f_grad=nothing,
                     objective_f_args=nothing,
                     constraint_names::Array{String,1}=["none"],
                     constraints_f=[(x) -> true],
                     constraints_f_grads=nothing,
                     constraints_f_args=nothing,
                     sample_points_budget=1000,
                     verbose=false, MAP_SIDE=800)::OptMapNLopt

    O = OptMapNLopt( variable_names,
                variable_domains,
                constraint_names,
                objective_f,
                objective_f_grad,
                objective_f_args,
                constraints_f,
                constraints_f_grads,
                constraints_f_args,
                sample_points_budget,
                verbose,
                zeros(1),
                zeros(1),
                zeros(1,1),
                zeros(1,1),
                zeros(1,1),
                zeros(RGB{Normed{UInt8, 8}}, MAP_SIDE, MAP_SIDE),
                zeros(RGB{Normed{UInt8, 8}}, MAP_SIDE, MAP_SIDE),
                nothing,
                MAP_SIDE)
    return O
end

"""
$(SIGNATURES)

Updates the optimization object as generates the heat map object

Arguments:
=========
O::OptMapNLopt:                             OptMapNLopt object
constraint_brightness::Float64=0.7:         setting the contrast factor in
                                            the brightness of the heatmap
                                            based on the constraint
                                            satisfiability
colormap=:viridi:                           colormap used for the heat map

Returns:
=======

"""
function create_map!(O::OptMapNLopt, constraint_brightness::Float64=0.7;
                     colormap=:viridis)
    # Generate the color map with the user-supplied colormap
    cmap = colorschemes[colormap]

    # Call the fit!() function to update the OptMapNLopt object
    fit!(O)

    # Extract the constraints and the arguments
    n_constraints = length(O.constraints_f)
    cargs         = O.constraints_f_args
    cgrads        = O.constraints_f_grads

    # Default behavior to generate the vectors of nothings
    if cargs === nothing
        args = [nothing for _ in 1:length(O.variable_domains)]
    end
    if cgrads === nothing
        cgrads = [nothing for _ in 1:length(O.variable_domains)]
    end

    # Extract the function arguments
    fargs = O.objective_f_args
    fgrad = O.objective_f_grad

    O.verbose && println("Creating OptMapNLopt")

    Z = zeros(Float64, O.MAP_SIDE, O.MAP_SIDE)
    A = ones(Float64, O.MAP_SIDE, O.MAP_SIDE)
    C = ones(Float64, O.MAP_SIDE^2, n_constraints)

    @inbounds for i in 1:O.MAP_SIDE
        @inbounds for j in 1:O.MAP_SIDE
            k = ((i - 1) * O.MAP_SIDE + j)

            @inbounds for c in 1:n_constraints
                if (cgrads[c] != nothing) && (cargs[c] != nothing)
                    C[k,c] = O.constraints_f[c](O.X_nd[:,k], cgrads[c],
                                                cargs[c]...)
                else
                    C[k,c] = O.constraints_f[c](O.X_nd[:,k])
                end
            end

            is_feasible_point = *(C[k,:]...) # product of all constraints
            if (fgrad != nothing) && (fargs != nothing)
                Z[i,j] = O.objective_f(O.X_nd[:,k], fgrad, fargs...)
            else
                Z[i,j] = O.objective_f(O.X_nd[:,k])
            end
            A[i,j] = Bool(is_feasible_point) ? 1.0 : constraint_brightness
        end
    end

    M = get(cmap, Z, (minimum(Z), maximum(Z)))
    M = RGB.(M .* A)
    M = convert(Array{RGB{N0f8}}, M)

    O.Z = Z
    O.M = M
    O.C = C
end

"""
$(SIGNATURES)

Draws the solution path from start to end vector projected onto the principal
axes from the ND space

Arguments:
=========

O::OptMapNLopt:
path_to_solution::AbstractArray{Float64,2}:
point_radius::Int=50:
line_width::Float64=3.0:

"""
function draw_solution!(O::OptMapNLopt,
        path_to_solution::AbstractArray{Float64,2}; point_radius::Int=50,
        line_width::Float64=3.0)
    O.verbose && println("Drawing OptMapNLopt")

    path_2d_tmp = MultivariateStats.transform(O.pca, path_to_solution)
    path_2d = zeros(Float64, 2, size(path_to_solution, 2))

    @inbounds for i in 1:size(path_to_solution, 2)
        path_2d[:,i] = coords_to_pixels(path_2d_tmp[:,i], O.min_2d, O.max_2d;
                                        maxgrid=O.MAP_SIDE)
    end

    tmp_map_name = tempname() * ".png"
    save(tmp_map_name, O.M)
    img = readpng(tmp_map_name)

    w = img.width
    h = img.height

    Drawing(w, h, tmp_map_name)
    placeimage(img, 0, 0, 1.0)
    setline(line_width)

    sethue("white")
    circle(path_2d[1,1], path_2d[2,1], point_radius, :stroke)
    sethue("red")
    circle(path_2d[1,end], path_2d[2,end], point_radius, :stroke)

    sethue(1.0, 0.0, 1.0)
    @inbounds for i in 2:size(path_2d, 2)
        p1 = Point((path_2d[:,i-1])...)
        p2 = Point((path_2d[:,i])...)

        circle(p1, line_width*2, :fill)
        circle(p2, line_width*2, :fill)
        line(p1, p2, :stroke)
    end

    finish()

    O.S = load(tmp_map_name)
    rm(tmp_map_name)
end

function save_base_img(O::OptMapNLopt, file_name::String)
    O.verbose && println("Saving image to $file_name")
    save(file_name, O.M)
end

function save_map(O::OptMapNLopt, file_prefix::String, path::String=".")
    map_name = "$(path)/$(file_prefix)_optmap.png"

    O.verbose && println("Saving map to $map_name")
    save(map_name, O.S)
end

function save_all(O::OptMapNLopt, file_prefix::String, path::String=".")
    save_map(O, file_prefix, path)

    X_csv = "$(path)/$(file_prefix)_optmap_X.csv"
    C_csv = "$(path)/$(file_prefix)_optmap_C.csv"
    Z_csv = "$(path)/$(file_prefix)_optmap_Z.csv"

    # O.verbose && println("Saving X_nd to $X_csv")
    CSV.write(X_csv, DataFrame(Matrix(O.X_nd'), :auto), header=O.variable_names)

    O.verbose && println("Saving C to $C_csv")
    CSV.write(C_csv, DataFrame(O.C, :auto), header=O.constraint_names)

    O.verbose && println("Saving Z to $Z_csv")
    CSV.write(Z_csv, DataFrame(O.Z, :auto), header=false)
end

function fit!(O::OptMapNLopt)
    O.verbose && println("Sampling nD points")
    X = sample_nd_points(O.variable_domains, O.sample_points_budget)

    O.verbose && println("Projecting nD -> 2D ($(size(X)))")
    X_2d, p = project_2d(X)

    O.verbose && println("Unprojecting 2D point grid")
    X_nd = inverse_project(X_2d, p, O.MAP_SIDE)

    O.pca = p
    O.X_nd = X_nd

    O.min_2d = [minimum(X_2d[1,:]), minimum(X_2d[2,:])]
    O.max_2d = [maximum(X_2d[1,:]), maximum(X_2d[2,:])]
end

function sample_nd_points(variables_domains::Array{Tuple{R,R},1} where R <:
        Number, point_budget::Int)::AbstractArray{Float64,2}
    n_vars = length(variables_domains)

    max_points = convert(Int, floor(point_budget^(1 / n_vars)))

    variable_ranges = []

    @inbounds for vd in variables_domains
        if typeof(vd[1]) == Float64
            points = max_points
        else
            points = vd[2] - vd[1] + 1
            points = min(points, max_points)
        end

        var_range = range(vd[1], vd[2], length=points)
        push!(variable_ranges, var_range)
    end

    tmp = vec(collect(product(variable_ranges...)))
    X = zeros(n_vars, length(tmp))

    @inbounds for j in 1:n_vars
        for i in 1:length(tmp)
            X[j,i] = tmp[i][j]
        end
    end

    return X
end

function project_2d(X::AbstractArray{Float64,2})
    p = fit(PCA, X; maxoutdim=2, method=:svd)
    X_2d = MultivariateStats.transform(p, X)

    return X_2d, p
end

function inverse_project(X_2d::AbstractArray{Float64,2}, p,
                         MAP_SIDE::Int64)::AbstractArray{Float64,2}
    tmp_grid = vec(collect(product(
                            range(minimum(X_2d[1,:]), maximum(X_2d[1,:]),
                                  length=MAP_SIDE), range(minimum(X_2d[2,:]),
                                                          maximum(X_2d[2,:]),
                                                          length=MAP_SIDE))))
    X_2d_grid = zeros(Float64, 2, MAP_SIDE^2)

    @inbounds for i in 1:size(tmp_grid)[1]
        X_2d_grid[1,i] = tmp_grid[i][1]
        X_2d_grid[2,i] = tmp_grid[i][2]
    end

    X_nd = reconstruct(p, X_2d_grid)
    return X_nd
end

function coords_to_pixels(point::AbstractArray{Float64,1},
        min_2d::Array{Float64,1}, max_2d::Array{Float64,1}; mingrid=1,
        maxgrid=100)
    p = (point - min_2d) ./ (max_2d - min_2d)
    p_trunc = floor.(p .* (maxgrid - mingrid))

    return p_trunc
end

end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
struct one_matching_output
    q::Float64
    mi::Vector{Int64}
end

struct all_matching_output
    q::Vector{Float64}
    mi::Vector{Int64}
    mj::Vector{Int64}
    medges::Int64
end

function intmatch(n::Int64,m::Int64,nedges::Int64,
                v1::Vector{Int64},v2::Vector{Int64},weight::Vector{Float64})

l1 = zeros(Float64,n)
l2 = zeros(Float64,n+m)
s = ones(Int64,n+m)
t = -1*ones(Int64,n+m)
offset = zeros(Int64,n)
deg = ones(Int64,n)
list = zeros(Int64,nedges+n)
index = zeros(Int64,nedges+n)
w = zeros(Float64,nedges+n)
match1 = -1*ones(Int64,n)
match2 = -1*ones(Int64,n+m)
tmod = zeros(Int64,m+n)
ntmod = 0

for i = 1:nedges
    deg[v1[i]] += 1
end

for i  = 2:n
    offset[i] = offset[i-1] + deg[i-1]
end
offset .+= 1

deg .= 0

for i  = 1:nedges
    list[offset[v1[i]] + deg[v1[i]]] = v2[i]
    w[offset[v1[i]] + deg[v1[i]]] = weight[i]
    index[offset[v1[i]] + deg[v1[i]]] = i
    deg[v1[i]] += 1
end

for i = 1:n
    list[offset[i]+deg[i]] = m + i
    w[offset[i]+deg[i]] = 0
    index[offset[i]+deg[i]] = -1
    deg[i] += 1
end

for i = 1:n
    for j = 0:deg[i]-1
        if w[offset[i]+j] > l1[i]
            l1[i] = w[offset[i]+j]
        end
    end
end

i = 1
while i<=n
    for j = 1:ntmod
        t[tmod[j]] = -1
    end
    ntmod = 0
    p = 1
    q = 1
    s[1] = i
    while p<=q
        if match1[i] >= 1
            break
        end
        k = s[p]
        for r = 0:deg[k]-1
            j = list[offset[k]+r]
            val1totest = w[offset[k]+r]
            val2totest = l1[k] + l2[j] - 1e-8
            if w[offset[k]+r] < val2totest
                continue
            end
            if t[j] < 0
                q = q + 1
                s[q] = match2[j]
                t[j] = k
                ntmod += 1
                tmod[ntmod] = j
                if match2[j] < 0
                    while j>=1
                        k = match2[j] = t[j]
                        p = match1[k]
                        match1[k] = j
                        j = p
                    end
                    break
                end
            end
        end
        p += 1
    end
    p -= 1

    if match1[i] < 0
        al = 1e20
        for j = 1:p
            t1 = s[j]
            for k = 0:deg[t1]-1
                t2 = list[offset[t1] + k]
                if t[t2] < 0 && l1[t1] + l2[t2] - w[offset[t1] + k] < al
                    al = l1[t1] + l2[t2] - w[offset[t1] + k]
                end
            end
        end
        for j = 1:p
            vtemp = s[j]
            vvtemp = l1[s[j]]
            l1[s[j]] -= al
            vvtemp = l1[s[j]]
        end
        for j = 1:ntmod
            l2[tmod[j]] += al
        end
    else
        i += 1
    end
end

ret = 0.0
for i = 1:n
    for j = 0:deg[i] - 1
        if list[offset[i] + j] == match1[i]
            ret += w[offset[i]+j]
        end
    end
end

mi = zeros(Int64,nedges)
for i = 1:n
    if match1[i]<=m
        for j = 0:deg[i]-1
            if list[offset[i] + j] == match1[i]
                mi[index[offset[i]+j]] = 1
            end
        end
    end
end

return one_matching_output(ret,mi)
end

function column_maxmatchsum(M::Int64,N::Int64,Qp::Vector{Int64},Qr::Vector{Int64},Qv::Vector{Float64},
                                m::Int64,n::Int64,nedges::Int64,li::Vector{Int64},lj::Vector{Int64})

minnm = n
if m<minnm
    minnm = m
end

q = zeros(Float64,N)
mi = zeros(Int64,Qp[N]-1)
mj = zeros(Int64,Qp[N]-1)

medges = 1

lwork1 = ones(Int64,m)
lind1 = -1*ones(Int64,m)
lwork2 = ones(Int64,n)
lind2 = -1*ones(Int64,n)

max_col_nonzeros = 0
for j = 1:N
    col_nonzeros = Qp[j+1] - Qp[j]
    if col_nonzeros >= max_col_nonzeros
        max_col_nonzeros = col_nonzeros
    end
end

se1 = zeros(Int64,max_col_nonzeros)
se2 = zeros(Int64,max_col_nonzeros)
sw = zeros(Float64,max_col_nonzeros)
sqi = zeros(Int64,max_col_nonzeros)
sqj = zeros(Int64,max_col_nonzeros)

for j = 1:N
    smalledges = 1
    nsmall1 = 1
    nsmall2 = 1

    for nzi = Qp[j]:Qp[j+1]-1
        i = Qr[nzi]
        v1 = li[i]
        v2 = lj[i]
        sv1 = -1
        sv2 = -1
        if lind1[v1] < 0
            sv1 = nsmall1
            lind1[v1] = sv1
            lwork1[sv1] = v1
            nsmall1 += 1
        else
            sv1 = lind1[v1]
        end
        if lind2[v2] < 0
            sv2 = nsmall2
            lind2[v2] = sv2
            lwork2[sv2] = v2
            nsmall2 += 1
        else
            sv2 = lind2[v2]
        end
        se1[smalledges] = sv1
        se2[smalledges] = sv2
        sw[smalledges] = Qv[nzi]
        sqi[smalledges] = i
        sqj[smalledges] = j
        smalledges += 1
    end
    nsmall1 -= 1
    nsmall2 -= 1
    smalledges -= 1

    if smalledges==0
        q[j] = 0
        continue
    end

    one_matching = intmatch(nsmall1,nsmall2,smalledges,se1,se2,sw)
    q[j] = one_matching.q
    smi = one_matching.mi

    for k = 1:smalledges
        if smi[k] > 0
            mi[medges] = sqi[k]
            mj[medges] = sqj[k]
            medges += 1
        end
    end

    for k = 1:nsmall1
        lind1[lwork1[k]] = -1
    end

    for k = 1:nsmall2
        lind2[lwork2[k]] = -1
    end

end

medges -= 1
return all_matching_output(q,mi,mj,medges)

end


text:
module ValidatedTransfer
using IntervalArithmetic, Compat, FastTransforms, JLD, FileIO

include("general.jl")
include("transforms.jl")
include("transferbounds.jl")

function filltransfer(J,K,transferfn_N,helperdata,aliaserror,entrybound,N_interp=maximum(J),T=BigInterval)
  S = Array{T}(length(J),length(K))
  for (k_ind,k) in enumerate(K)
    B = Array{T}(N_interp)
    for i in eachindex(B)
        B[i] = transferfn_N(i,k,helperdata)
    end
    S[:,k_ind] = chebyshevtransform(B)[J]
    for j in J
      ae = aliaserror(j,k,N_interp)
      S[j,k_ind] += erroradjustment(Interval(ae))
      S[j,k_ind] = S[j,k_ind] ∩ erroradjustment(entrybound(j,k))
    end
  end

  S
end

function maketransfer(N,transferfn,helperdata,aliaserror,entrybound,workerprocs=workers(),N_interp=N,T=BigInterval)
  nws = length(workerprocs)
  enum_workers = enumerate(workerprocs)

  M = Array(T,N,N)
  @sync begin
    for (pn,pid) in enum_workers
      @async M[:,pn:nws:N] = remotecall_fetch(filltransfer,pid,1:N,pn:nws:N,transferfn,helperdata,aliaserror,entrybound,N_interp,T)
    end
  end
  M
end

function makesolutioninv!(M::Array)
  # M = fetch(M_r)
  for k in 1:size(M,2)
    for j = 1:size(M,1)
      M[j,k] *= -1
    end

    mod(k,2) == 1 && (k== 1 ? (M[1,k] += 1) :
        (M[1,k] -= 1/convert(eltype(M),(k-1)^2-1)) )
    M[k,k] += 1
  end
  M
end

include("inversion.jl")

function getsolution!(M,workerprocs=workers();verbose=true)
  makesolutioninv!(M)
  parallel_invert(M,workerprocs;verbose=verbose)
end

function getsolution(M,workerprocs=workers();verbose=true)
  K=makesolutioninv!(deepcopy(M))
  parallel_invert(K,workerprocs;verbose=verbose)
end

include("files.jl")

export BigInterval, erroradjustment
export ChebyshevZeroOneBV, AnalyticEntryBounds, FirstOrderLYBounds
export logSEnorm, logEnorm, Snorm, aliaserror, entrybound
export transfer_basisfun, transfer_basisfun_multiel
export maketransfer, getsolution
export savezip, loadzip

end # module


text:
"""
    SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <: AbstractBridgeOptimizer

The `SingleBridgeOptimizer` bridges any objective functions supported by the
bridge `BT`. This is in contrast with the [`MathOptInterface.Bridges.LazyBridgeOptimizer`](@ref)
which only bridges the objective functions that are unsupported by the internal model,
even if they are supported by one of its bridges.
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge,OT<:MOI.ModelLike} <:
               MOIB.AbstractBridgeOptimizer
    model::OT
    map::Map # `MOI.ObjectiveFunction` -> objective bridge
end
function SingleBridgeOptimizer{BT}(model::OT) where {BT,OT<:MOI.ModelLike}
    return SingleBridgeOptimizer{BT,OT}(model, Map())
end

function bridges(bridge::MOI.Bridges.AbstractBridgeOptimizer)
    return EmptyMap()
end
bridges(bridge::SingleBridgeOptimizer) = bridge.map

MOIB.supports_constraint_bridges(::SingleBridgeOptimizer) = false
function MOIB.is_bridged(::SingleBridgeOptimizer, ::Type{<:MOI.AbstractSet})
    return false
end
function MOIB.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end
function MOIB.supports_bridging_objective_function(
    ::SingleBridgeOptimizer{BT},
    F::Type{<:MOI.AbstractScalarFunction},
) where {BT}
    return supports_objective_function(BT, F)
end
function MOIB.is_bridged(
    b::SingleBridgeOptimizer,
    F::Type{<:MOI.AbstractScalarFunction},
)
    return MOIB.supports_bridging_objective_function(b, F)
end
function MOIB.bridge_type(
    ::SingleBridgeOptimizer{BT},
    ::Type{<:MOI.AbstractScalarFunction},
) where {BT}
    return BT
end


text:
const alphavantage_api = "https://www.alphavantage.co/"

struct AVClient
    key::String
    entry::String
end

AVClient(; key = "", entry = alphavantage_api) = AVClient(key, entry)

const GLOBAL = Ref(AVClient(key = get(ENV, "ALPHA_VANTAGE_API_KEY", "")))

function global_key!(key)
    GLOBAL[] = AVClient(key, GLOBAL[].entry)
end

function global_entry!(entry)
    GLOBAL[] = AVClient(GLOBAL[].key, entry)
end

function key(client::AVClient)
    if isempty(client.key)
        @warn "No API key found"
    end

    return client.key
end

entry(client::AVClient) = client.entry


text:
if(!isdefined(:Calculus))           using Calculus        end
if(!isdefined(:ForwardDiff))        using ForwardDiff       end

# Manually add code to autodiff mod2pi
function Calculus.differentiate(::Calculus.SymbolParameter{:mod2pi}, args, wrt)
   x = args[1]
   xp = Calculus.differentiate(x, wrt)
   if xp != 0
       return Calculus.@sexpr(xp)
   else
       return 0
   end
end

function Base.mod2pi(n::ForwardDiff.Dual)
   return ForwardDiff.Dual(mod2pi(ForwardDiff.value(n)), ForwardDiff.partials(n))
end


text:
using Mocking

send_sshpublic_key_1_patch = @patch function send_sshpublic_key(a...; b...)
    return Dict{String, Any}("Success" => true, "RequestId" => "abcd1234-abcd-1234-abcd-1234abcd1234")
end


text:
# detect if using SLURM
const IN_SLURM = "SLURM_JOBID" in keys(ENV)

IN_SLURM && using ClusterManagers

using ShaleDrillingModel
using Test
using StatsFuns
using JLD

# load in prices, costs, and transition between high/low volatility regimes
# see eample/make_big_transition.jl
jldpath = joinpath(Pkg.dir("ShaleDrillingModel"), "data/price-cost-transitions.jld")
@load jldpath pspace cspace Πpcr Πp1

# -------- state space ----------

# set of observable types
royalty_rates = [0.125, 1./6., 0.1875, 0.20, 0.225, 0.25]
geology_types = 1:10

# number of points in unobserved heterogeneity grid
nψ, nv =  5, 5

# grids for prices & unobserved heterogeneity
zspace, ψspace, vspace = (pspace,), range(-5.5, stop=5.5, length=nψ),  range(-3.0, stop=3.0, length=nv)

# Structure of time / number of wells
# drill 0 to 2 wells/period
# can drill up to 4 wells total
# 10+1 periods before lease expires
wp = well_problem(2, 4, 10)

# discount rate
β = (1.02 / 1.125) ^ (1./12.)  # real discount rate

# ------------- set up primitives & tmp vars -------------

# initial parameters
θt = [1.0, -2.6, -2.4, 1.2]  # flow payoff parameters
σv = 0.5                     # unobserved heterogeneity

# primitives holds (almost) all information to solve VF except observable type
# flow payoffs + gradients are u_add, du_add, and duσ_add
prim = dcdp_primitives(u_add, du_add, duσ_add, β, wp, zspace, Πp1, ψspace, vspace, 1)

# temporary variables
tmpv = dcdp_tmpvars(length(θt), prim)

# Emax + gradients
evs = dcdp_Emax(θt, prim)

# check sizes of primitives
ShaleDrillingModel.check_size(θt, prim, evs)

# check flow gradients w/ finite differences
check_flowgrad(θt, σv, prim, 0.2)

# solve problem & check jacobian of Emax w/ finite difference
check_EVjac(evs, tmpv, prim, θt, σv, 0.2)

# ----------- how we do this in parallel ------------

# set up workers
pids = IN_SLURM ? addprocs_slurm(parse(Int, ENV["SLURM_CPUS_PER_TASK"])) : addprocs()

# tell workers about the pkg
@everywhere @show pwd()
@everywhere using ShaleDrillingModel

# initialize struct of shared arrays
sev = SharedEV(pids, vcat(θt, σv), prim, royalty_rates, 1:1)

# send primitives to workers
@eval @everywhere begin
    set_g_dcdp_primitives($prim)
    set_g_dcdp_tmpvars($tmpv)
    set_g_SharedEV($sev)
end

# parallel solve
parallel_solve_vf_all!(sev, vcat(θt,σv), Val{true})

using Plots
gr()

nSexpl = ShaleDrillingModel._nSexp(prim)  # number of states in initial period
nS = ShaleDrillingModel._nS(prim)         # number of endogenous states overall

plot(prim.wp.τmax:-1:0, sev.EV[15, 3, 1:nSexpl, :, 1], xaxis=(:flip,), label=string.(round.(royalty_rates,3)), xlabel="time remaining", ylabel="Emax(V)")
plot(exp.(pspace), sev.EV[:, 3, nSexpl+1:nS, 4, 1], xlabel="Price", ylabel="Emax(V)")


text:
module ScatteringTransform

using ChainRules
using CUDA
using FFTW
using Wavelets, ContinuousWavelets
using Zygote, Flux, LinearAlgebra, AbstractFFTs
using Flux
using FourierFilterFlux
using Adapt
using RecipesBase
using Base: tail
using ChainRulesCore
using Plots # who cares about weight really?
using Statistics
using Dates

import Adapt: adapt
import ChainRules: rrule
import Zygote: has_chain_rrule, rrule
import Wavelets: eltypes

include("shared.jl")
export scatteringTransform, stFlux, stParallel, eltypes, depth
include("pathLocs.jl")
export pathLocs
include("scattered.jl")
export Scattered, ScatteredFull, ScatteredOut, nonZeroPaths, cat
include("pool.jl")
export RationPool, nPoolDims, outdims, poolSize
include("transform.jl")
export normalize, cu

function scatteringTransform(inputSize, m, backend::UnionAll; kwargs...)
    backend(inputSize, m; kwargs...)
end

scatteringTransform(inputSize, m; kwargs...) =
    scatteringTransform(inputSize, m, stFlux; kwargs...)

include("utilities.jl")
export getWavelets, flatten, roll, importantCoords, batchOff, getParameters, getMeanFreq
export roll, wrap, flatten
include("adjoints.jl")
include("interpretationTools.jl")
export ∇st, plotFirstLayer1D, gifFirstLayer, plotSecondLayer, addNextPath, jointPlot
end # end Module


text:
"""
    setup_sample_jobs(genome_path, genome_index_path, gff3_path, sample_set_length, sample_window_min, sample_window_max, perigenic_pad; deterministic)

Return sample job channels for 'execute_sample_jobs', given 'genome_path' to a properly formatted FASTA genome, 'genome_index_path' to the associated .FAI, 'gff3_path' to the genome's GFF3 feature database, a total number of bases to sample, 'sample_set_length', as well as 'sample_window_min' and 'sample_window_max' lengths for individual samples. 'perigenic_pad' specifies the number of bases up- and down-stream of exonic sequences to be considered 'perigenic'; if 0, the perigenic partition will be entirely intronic. Optional boolean 'deterministic' may be set to 'true' to always return '+' strand sequences from unstranded samples; otherwise unstranded samples may be returned with either orientation.
"""

#function to partition genome and set up Distributed RemoteChannels so partitions can be sampled simultaneously
function setup_sample_jobs(genome_path::String, genome_index_path::String, gff3_path::String, sample_set_length::Integer, sample_window_min::Integer, sample_window_max::Integer, perigenic_pad::Integer; deterministic::Bool=false)
    #argument checking
    !ispath(genome_path) && throw(ArgumentError("Bad genome path!"))
    !ispath(genome_index_path) && throw(ArgumentError("Bad genome index path!"))
    !ispath(gff3_path) && throw(ArgumentError("Bad gff3 path!"))
    sample_set_length < 1 && throw(ArgumentError("Sample set length must be a positive integer!"))
    sample_window_min < 1 || sample_window_max < 1 && throw(ArgumentError("Sample window minimum and maximum bounds must be positive integers!"))
    sample_window_min >= sample_window_max && throw(ArgumentError("Sample window minimum size must be smaller than maximum size"))
    perigenic_pad < 0 && throw(ArgumentError("Perigenic pad must be 0 or positive!"))

    coordinate_partitions = partition_genome_coordinates(gff3_path, perigenic_pad)
    sample_sets = DataFrame[]
    input_sample_jobs = RemoteChannel(()->Channel{Tuple}(length(coordinate_partitions))) #channel to hold sampling jobs
    completed_sample_jobs = RemoteChannel(()->Channel{Tuple}(length(coordinate_partitions))) #channel to hold completed sample dfs
    for (partitionid, partition) in coordinate_partitions
        add_metacoordinates!(partition)
        put!(input_sample_jobs, (genome_path, genome_index_path, partition, partitionid, sample_set_length, sample_window_min, sample_window_max, deterministic))
    end
    progress_channel = RemoteChannel(()->Channel{Tuple}(20))
    return (input_sample_jobs, completed_sample_jobs, progress_channel, sample_set_length)
end

"""
    execute_sample_jobs!(channels, worker_pool)

Return sample dataframes for each genomic partition, given appropriately set up job channels (from setup_sample_jobs!) and a vector of worker ids ('worker_pool', [1] for master).
"""

function execute_sample_jobs(channels::Tuple{RemoteChannel,RemoteChannel,RemoteChannel,Integer}, worker_pool::Vector{Int64}; partitions::Integer=3)
    input_sample_channel, completed_sample_channel, progress_channel, sample_set_length = channels

    #argument checking
    length(worker_pool) < 1 && throw(ArgumentError("Worker pool must contain one or more worker IDs!"))
    !isready(input_sample_channel) && throw(ArgumentError("Input sample channel not ready, likely channel set from setup_sample_jobs passed incorrectly"))
    sample_set_length < 1 && throw(ArgumentError("Sample set length must be a positive integer, likely channel set from setup_sample_jobs passed incorrectly"))

    #send sampling jobs to workers
    if isready(input_sample_channel) > 0
        @info "Sampling.."
        #WORKERS SAMPLE
        for (n, worker) in enumerate(worker_pool)
            if n <= partitions #no more workers than partitions to be used
                remote_do(get_sample_set, worker, input_sample_channel, completed_sample_channel, progress_channel)
            end
        end
    else
        @error "No sampling jobs!"
    end

    #progress meters for sampling
    sampling_meters=Dict{String, Progress}()
    overall_sampling_meter=Progress(partitions,"Overall sampling progress:")
    completed_counter = 0
    ProgressMeter.update!(overall_sampling_meter, completed_counter)
    sampling_offset = ones(Bool, partitions)

    #collect progress updates while waiting on completion of sampling jobs
    while completed_counter < partitions
        wait(progress_channel)
        partition_id, progress = take!(progress_channel)
        if haskey(sampling_meters, partition_id)
            ProgressMeter.update!(sampling_meters[partition_id], progress)
        else
            offset = findfirst(sampling_offset)[1]
            sampling_meters[partition_id] = Progress(sample_set_length, "Sampling partition $partition_id:", offset)
            ProgressMeter.update!(sampling_meters[partition_id], progress)
            sampling_offset[offset] = false
        end
        if progress == sample_set_length
            completed_counter += 1
            ProgressMeter.update!(overall_sampling_meter, completed_counter)
        end
    end

    #collect sample dfs by partition id when ready
    sample_record_dfs = Dict{String,DataFrame}()
    collected_counter = 0
    while collected_counter < partitions
        wait(completed_sample_channel)
        partition_id, sample_df = take!(completed_sample_channel)
        sample_record_dfs[partition_id] = sample_df
        collected_counter += 1
        @info "Partition $partition_id completed sampling..."
    end

    return sample_record_dfs
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# simulation replications

# set simulation replications, one rep per call set
# mutates: sim, callSets
function setSimReps!(sim::Simulation, callSets::Vector{Vector{Call}})
	sim.reps = []
	for calls in callSets
		initCalls!(sim, calls)
		rep = Simulation()
		rep.calls = calls
		rep.numCalls = length(calls)
		push!(sim.reps, rep)
	end
	sim.numReps = length(sim.reps)
	if isdefined(sim, :backup) sim.backup.numReps = sim.numReps end

	# have not yet handled writing output files for each replication
	sim.writeOutput = false
	if isdefined(sim, :backup) sim.backup.writeOutput = false end
end

# run single simulation replication on `sim`
# need this function in case rep.isRunnable == false
# mutates: sim, rep
function simulateRep!(sim::Simulation, rep::Simulation)
	sim.calls = Call[] # so calls will not be reset
	reset!(sim)
	setSimCalls!(sim, rep.calls)
	simulate!(sim)
	for fname in (:startTime, :time, :endTime, :stats, :used, :complete)
		setfield!(rep, fname, getfield(sim, fname))
	end
end

# run simulation replications on `sim`
# need this function in case rep.isRunnable == false
# mutates: sim, reps
function simulateReps!(sim::Simulation, reps::Vector{Simulation} = sim.reps; doPrint::Bool = false)
	numReps = length(reps)
	for (i,rep) in enumerate(reps)
		doPrint && print("\rSimulating replication $i of $numReps.")
		simulateRep!(sim, rep)
	end
	doPrint && println()
end

# make sim replications runnable, e.g. in order to run simulate!(sim.reps[1])
# for experiments, any changes made to `sim` may require the same change to be made to each replication
# mutates: sim.reps
function makeRepsRunnable!(sim::Simulation)
	@assert(sim.initialised && !sim.used)
	for rep in sim.reps
		fnamesDontCopy = (:backup, :net, :travel, :grid, :resim, :calls, :demand, :demandCoverage, :reps)
		for fname in setdiff(fieldnames(Simulation), fnamesDontCopy)
			setfield!(rep, fname, deepcopy(getfield(sim, fname)))
		end

		# have some fields in rep point to same data as sim
		for fname in (:net, :grid, :demandCoverage)
			setfield!(rep, fname, getfield(sim, fname))
		end

		# travel and demand types have some fields that change value while simulating
		# in particular: sim.travel.recentSetsStartTimesIndex (Int) and sim.demand.recentSetsStartTimesIndex (Int)
		for fname in (:travel, :demand)
			dst = getfield(rep, fname)
			src = getfield(sim, fname)
			for fname2 in fieldnames(typeof(dst))
				setfield!(dst, fname2, getfield(src, fname2))
			end
		end

		backup!(rep)
		setSimCalls!(rep, rep.calls) # requires rep.backup to be defined

		rep.isRunnable = rep.backup.isRunnable = true
	end
end

function simulateRep!(rep::Simulation)
	@assert(rep.isRunnable)
	simulate!(rep)
end

# run simulation replications
# requires replications to be 'runnable', see function makeRepsRunnable!()
# if parallel=true, will use multi-threading
# mutates: reps
function simulateReps!(reps::Vector{Simulation}; parallel::Bool = false, numThreads::Int = maxNumThreads, doPrint::Bool = false)
	@assert(allunique(reps))
	@assert(all(rep -> rep.isRunnable, reps))
	if parallel
		doPrint = false # IO not thread safe
		runParallel!(simulateRep!, reps...; numThreads = numThreads)
	else
		numReps = length(reps)
		for (i,rep) in enumerate(reps)
			doPrint && print("\rSimulating replication $i of $numReps.")
			simulateRep!(rep)
		end
		doPrint && println()
	end
end

function resetRep!(rep::Simulation)
	@assert(rep.isRunnable)
	reset!(rep)
end

function resetReps!(reps::Vector{Simulation}; parallel::Bool = false, numThreads::Int = maxNumThreads)
	@assert(allunique(reps))
	@assert(all(rep -> rep.isRunnable, reps))
	if parallel
		runParallel!(resetRep!, reps...; numThreads = numThreads)
	else
		for rep in reps resetRep!(rep) end
	end
end


text:
"""
Plots history of sampled points or init_points
"""
function visualize_search(hist, color = "blue")
    n = length(hist)
    xx = zeros(1, n)
    yy = zeros(1, n)

    for i = 1:n
        xx[i] = hist[i][1]
    end
    for i = 1:n
        yy[i] = hist[i][2]
    end
    Plots.scatter!(xx, yy, color = color, legend = false)
end

function visualize_search2(hist, color = "red")
    n = size(hist, 2)
    xx = zeros(1, n)
    yy = zeros(1, n)

    for i = 1:n
        xx[i] = hist[i][1]
    end
    for i = 1:n
        yy[i] = hist[i][2]
    end
    Plots.scatter!(xx, yy, color = color, legend = false)
end
visualize_search(init_hist, "red")
visualize_search(x_hist)
#labels = ["1", "2", "3"]

#for i = 1:25
#    init_hist[i] = init_hist[i][2:end]
#end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module QCCDevControl_Types
export GateZone, Junction, AuxZone, Qubit, LoadingZone, QubitStatus, typesSizes, JunctionType
export Time_t, QCCDevControl

using LightGraphs
using ..QCCDevDes_Types

#Types of zones
const zoneTypes = Set([:junction,:loadingZone, :auxZone, :gateZone])

# Possible Qubits Status (not used)
const QubitStatus = Set([:inLoadingZone, :inGateZone])

# Supported junction types with corresponding sizes
const JunctionType = Set([:T, :Y, :X ])
const typesSizes = Dict(:T => 3, :Y => 3, :X => 4)

"""
Type for time inside the qdev, in [change if necessary]   10^{-10}
seconds, i.e., ns/10.  All times are ≥0; negative value of expressions
of this type are errors (and may carry local error information).
"""
const Time_t = Int64

"""
Struct for junction.
id: Junction ID.
type: Type of the junction. Each type may define how the junction works differently
ends: Array of its connections Id's
Throws ArgumentError if junction type doesn't match with number of ends.
"""
struct Junction
    id::Symbol
    type::Symbol
    zoneType::Symbol
    ends::Array{Symbol}
    function Junction(id::Symbol, type::Symbol, ends::Array{Symbol})
        type in JunctionType || throw(ArgumentError("Junction type $type not supported"))
        if length(ends) != typesSizes[type]
            throw(ArgumentError("Junction with ID $id of type $type has $(length(ends)) ends." *
            " It should have $(typesSizes[type]) ends."))
        end
        return new(id, type, :junction, ends)
    end
end

"""
Struct for the qubits.
id: qubit ID
status: current qubit status
    - moving
    - resting
    - waitingDecongestion
    - gateApplied
position: current qubit position
destination: qubit destination, it could not have any
"""
mutable struct Qubit
    id::Int
    status::Symbol
    position::Symbol
    destination::Union{Nothing,Symbol}
    Qubit(id::Int, position::Symbol) = new(id, :inLoadingZone, position, nothing)
end

"""
Struct for the auxiliary zones.
id: trap identifier
capacity: maximum qubits in the trap
chain: Ordered Qubits in the trap (from end0 to end1)
end0 & end1: Trap endings
"""
mutable struct LoadingZone
    id::Symbol
    zoneType::Symbol
    end0::Union{Symbol, Nothing}
    end1::Union{Symbol, Nothing}
    hole::Union{Int, Nothing}
    LoadingZone(id, end0, end1) = end0 == end1 && (!isnothing(end0) || !isnothing(end1)) ?
    throw(ArgumentError("In loading zone $id : \"end0\" and \"end1\" must be different")) :
    new(id, :loadingZone, end0, end1, nothing)
end

"""
Struct for the auxiliary zones.
id: trap identifier
capacity: maximum qubits in the trap
chain: Ordered Qubits in the trap (from end0 to end1)
end0 & end1: Trap endings
"""
struct AuxZone
    id::Symbol
    zoneType::Symbol
    capacity::Int64
    end0::Union{Symbol, Nothing}
    end1::Union{Symbol, Nothing}
    chain::Array{Array{Int64,1},1}
    AuxZone(id, capacity, end0, end1) = end0 == end1 && (!isnothing(end0) || !isnothing(end1)) ?
    throw(ArgumentError("In aux zone $id : \"end0\" and \"end1\" must be different")) :
    new(id, :auxZone, capacity, end0, end1, [])
end

"""
Struct for the gate zones.
id: trap identifier
capacity: maximum qubits in the trap
chain: Ordered Qubits in the trap (from end0 to end1)
end0 & end1: Trap endings
"""
struct GateZone
    id::Symbol
    zoneType::Symbol
    capacity::Int64
    end0::Union{Symbol, Nothing}
    end1::Union{Symbol, Nothing}
    chain::Array{Array{Int64,1},1}
    GateZone(id, capacity, end0, end1) = end0 == end1 && (!isnothing(end0) || !isnothing(end1)) ?
    throw(ArgumentError("In gate zone $id : \"end0\" and \"end1\" must be different")) :
    new(id, :gateZone, capacity, end0, end1, [])
end

mutable struct QCCDevControl
    dev            ::QCCDevDescription

    simulate       ::Symbol                   # one of `:No`, `:PureStates`, `:MixedStates`
    qnoise_esimate ::Bool                     # whether estimation of noise takes place

    t_now          ::Time_t

    qubits      ::Dict{Int,Qubit}
    gateZones      ::Dict{Symbol,GateZone}
    junctions      ::Dict{Symbol,Junction}
    auxZones       ::Dict{Symbol,AuxZone}
    loadingZones   ::Dict{Symbol,LoadingZone}

    # graph          ::SimpleGraph{Int64}

    QCCDevControl(dev, simulate, qnoise_estimate, gateZones, junctions, auxZones, loadingZones) =
            new(dev, simulate, qnoise_estimate, 0, Dict{String,Qubit}(),
            gateZones, junctions, auxZones, loadingZones)

end

end

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
using GLVisualize, GeometryTypes, FileIO
using GLAbstraction, Colors, Reactive

if !isdefined(:runtests)
	window = glscreen()
	timesignal = loop(linspace(0f0,1f0,360))
end

cube = HyperRectangle(Vec3f0(0), Vec3f0(0.05))
n = 20
const wx,wy,wz = widths(cube)

mesh = GLNormalMesh(cube)

timepi = const_lift(*, timesignal, 2f0*pi)
function position(t, x, y)
    pos = Point3f0(x*(sqrt(wx^2+wy^2)), -y*wy, y*wz)
    dir = Point3f0(0, wy, wz)
    pos = pos + sin(t)*dir
end
position_signal = map(timepi) do t
    vec(Point3f0[position(t,x,y) for x=1:n, y=1:n])
end

rotation = map(timepi) do t
    vec(Vec3f0[Vec3f0(cos(t+(x/7)),sin(t+(y/7)), 1) for x=1:20, y=1:20])
end

cubes = visualize(
    (mesh, position_signal),
    rotation=rotation,
    color_map=GLVisualize.default(Vector{RGBA}),
    color_norm=Vec2f0(1,1.8)
    # intensity that will define the color sampled from color_map will fallback
    # to the length of the rotation vector.
    # you could also supply it via intensity = Vector{Float32}
)

# we create our own camera to better adjust to what we want to see.
camera = PerspectiveCamera(
    Signal(Vec3f0(0)), # theta (rotate by x around cam xyz axis)
    Signal(Vec3f0(0)), # translation (translate by translation in the direction of the cam xyz axis)
    Signal(Vec3f0(wx*n+4wx,wy*n,wz*n)/2), # lookat. We want to look at the middle of the cubes
    Signal(Vec3f0(((wx*n+4wx)/2),1.2*wy*n,(wz*n)/2)), # camera position. We want to be on the same height, but further away in y
    Signal(Vec3f0(0,0,1)), #upvector
    window.area, # window area

    Signal(41f0), # Field of View
    Signal(1f0),  # Min distance (clip distance)
    Signal(100f0), # Max distance (clip distance)
    Signal(GLAbstraction.ORTHOGRAPHIC)
)
_view(cubes, window, camera=camera)

if !isdefined(:runtests)
	renderloop(window)
end


text:
function update!(events::Events{T},
                 event::Event{T}) where T <: EpidemicModel
  events[event.new_state][event.individual] = event.time
  return events
end

function update!(states::Vector{DiseaseState},
                 event::Event{T}) where T <: EpidemicModel
  states[event.individual] = event.new_state
  return states
end

function update!(net::TransmissionNetwork,
                 xfer::ExogenousTransmission)
  net.external[xfer.individual] = true
  return net
end

function update!(net::TransmissionNetwork,
                 xfer::EndogenousTransmission)
  net.internal[xfer.source, xfer.individual] = true
  return net
end

function update!(net::TransmissionNetwork,
                 xfer::NoTransmission)
  return net
end

function update!(tr::TransmissionRates,
                 event::Event{T},
                 states::Vector{DiseaseState},
                 pop::Population,
                 rf::RiskFunctions{T},
                 rp::RiskParameters{T};
                 xzero::Bool=false,
                 transmission_rate_cache::Union{Nothing, TransmissionRateCache}=nothing) where T <: EpidemicModel
  id = event.individual
  if _new_transmission(event) && xzero
    tr.external[id] = 0.0
    tr.internal[:, id] .= 0.0
  end
  if event.new_state == State_I
    if transmission_rate_cache == nothing
      transmissibility = rf.transmissibility(rp.transmissibility, pop, id)
      @inbounds @simd for i in findall(states .== Ref(State_S))
        tr.internal[id, i] = rf.susceptibility(rp.susceptibility, pop, i) *
                                    rf.infectivity(rp.infectivity, pop, i, id) *
                                    transmissibility
      end
    else
      @inbounds @simd for i in findall(states .== Ref(State_S))
        if transmission_rate_cache.internal[id, i] == nothing
          tr.internal[id, i] = rf.susceptibility(rp.susceptibility, pop, i) *
                               rf.infectivity(rp.infectivity, pop, i, id) *
                               rf.transmissibility(rp.transmissibility, pop, id)
          transmission_rate_cache.internal[id, i] = tr.internal[id, i]
        else
          tr.internal[id, i] = transmission_rate_cache.internal[id, i]
        end
      end
    end
  elseif event.new_state == State_R
    tr.internal[id, states .== Ref(State_S)] .= 0.0
  end
  return tr
end

function update!(rates::EventRates{T},
                 tr::TransmissionRates,
                 event::Event{T},
                 states::Vector{DiseaseState},
                 pop::Population,
                 rf::RiskFunctions{T},
                 rp::RiskParameters{T}) where T <: EpidemicModel
  id = event.individual
  if event.new_state == State_E
    rates.exposure[id] = 0.0
    rates.infection[id] = rf.latency(rp.latency, pop, id)
    @logmsg LogLevel(-5000) "Exposure rate total for i = $id updated"  λ = rates.exposure[id]
    @logmsg LogLevel(-5000) "Infection rate for i = $id updated" λ = rates.infection[id]
  elseif event.new_state == State_I
    rates.infection[id] = 0.0
    @logmsg LogLevel(-5000) "Infection rate for i = $id updated" λ = rates.infection[id]
    if T in [SEIR; SIR]
      rates.removal[id] = rf.removal(rp.removal, pop, id)
      @logmsg LogLevel(-5000) "Removal rate for i = $id updated" λ = rates.removal[id]
    end
    if T in [SEIR; SEI]
      @simd for i in findall(states .== Ref(State_S))
        # This assumes `TransmissionRates` already updated!
        rates.exposure[i] += tr.internal[id, i]
        @logmsg LogLevel(-5000) "Exposure rate total for i = $id updated"  λ = rates.exposure[id]
      end
    elseif T in [SIR; SI]
      @simd for i in findall(states .== Ref(State_S))
        # This assumes `TransmissionRates` already updated!
        rates.infection[i] += tr.internal[id, i]
        @logmsg LogLevel(-5000) "Infection rate total for i = $id updated"  λ = rates.infection[id]
      end
    end
  elseif event.new_state == State_R
    rates.removal[id] = 0.0
    @logmsg LogLevel(-5000) "Removal rate for i = $id updated" λ = rates.removal[id]
    if T == SEIR
      @simd for i in findall(states .== Ref(State_S))
        # This assumes `TransmissionRates` already updated!
        rates.exposure[i] = sum(tr, i)
        @logmsg LogLevel(-5000) "Exposure rate total for i = $i updated"  λ = rates.exposure[i]
      end
    elseif T == SIR
      @simd for i in findall(states .== Ref(State_S))
        # This assumes `TransmissionRates` already updated!
        rates.infection[i] = sum(tr, i)
        @logmsg LogLevel(-5000) "Infection rate for i = $i updated" λ = rates.infection[i]
      end
    end
  end
  return rates
end

function update!(tr::TransmissionRates,
                 rates::EventRates{T},
                 event::Event{T},
                 states::Vector{DiseaseState},
                 pop::Population,
                 rf::RiskFunctions{T},
                 rp::RiskParameters{T};
                 transmission_rate_cache::Union{Nothing, TransmissionRateCache}=nothing) where T <: EpidemicModel
  update!(tr, event, states, pop, rf, rp, transmission_rate_cache=transmission_rate_cache)
  update!(rates, tr, event, states, pop, rf, rp)
  return tr, rates
end