text:
# Test basic extensions
@testset "Basic Extensions" begin
    # initialize the model, references, and other info
    m = InfiniteModel()
    cref1 = FiniteConstraintRef(m, 1, ScalarShape())
    cref2 = InfiniteConstraintRef(m, 2, ScalarShape())
    con = BoundedScalarConstraint(zero(AffExpr), MOI.EqualTo(0.0),
                                  ParameterBounds(), ParameterBounds())
    # owner_model
    @testset "JuMP.owner_model" begin
      @test owner_model(cref1) == m
    end
    # index
    @testset "JuMP.index" begin
      @test JuMP.index(cref1) == 1
    end
    # is_valid
    @testset "JuMP.is_valid" begin
      @test !is_valid(m, cref1)
      m.constrs[1] = con
      @test is_valid(m, cref1)
    end
    # constraint_object
    @testset "JuMP.constraint_object" begin
      @test constraint_object(cref1).func == con.func
      @test constraint_object(cref1).set == con.set
      @test constraint_object(cref1).bounds == con.bounds
      @test constraint_object(cref1).orig_bounds == con.orig_bounds
    end
    # test Base.:(==) of constraint references
    @testset "Base.:(==) References" begin
        @test !(cref1 == cref2)
        @test cref1 == FiniteConstraintRef(m, 1, ScalarShape())
        @test cref1 != FiniteConstraintRef(m, 2, ScalarShape())
        @test cref1 != FiniteConstraintRef(InfiniteModel(), 1, ScalarShape())
    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) == GeneralConstraintRef
    end
    # test shape
    @testset "shape" begin
        @test JuMP.shape(con) == ScalarShape()
    end
    # test jump_function
    @testset "JuMP.jump_function" begin
        @test jump_function(con) == zero(AffExpr)
    end
    # test moi_set
    @testset "JuMP.moi_set" begin
        @test moi_set(con) == MOI.EqualTo(0.0)
    end
end

# Test name methods
@testset "Name" 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)
    data = DiscreteMeasureData(par, [1], [1])
    meas = measure(inf + par - x, data)
    cref = FiniteConstraintRef(m, 1, ScalarShape())
    m.constr_to_name[1] = "test"
    # JuMP.name
    @testset "JuMP.name" begin
        @test name(cref) == "test"
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        @test isa(set_name(cref, "new"), Nothing)
        @test name(cref) == "new"
    end
    # _make_constraint_ref
    @testset "_make_constraint_ref" begin
        # prepare for infinite constraint
        m.constrs[1] = ScalarConstraint(inf + par, MOI.LessThan(1.0))
        @test InfiniteOpt._make_constraint_ref(m, 1) == InfiniteConstraintRef(m,
                                                        1, ScalarShape())
        # prepare for measure constraint
        m.constrs[1] = ScalarConstraint(pt + meas, MOI.LessThan(1.0))
        @test InfiniteOpt._make_constraint_ref(m, 1) == MeasureConstraintRef(m,
                                                        1, ScalarShape())
        # prepare for finite constraint
        m.constrs[1] = ScalarConstraint(pt + x, MOI.LessThan(1.0))
        @test InfiniteOpt._make_constraint_ref(m, 1) == FiniteConstraintRef(m,
                                                        1, ScalarShape())
    end
    # constraint_by_name
    @testset "JuMP.constraint_by_name" begin
        # test normal
        @test constraint_by_name(m, "new") == cref
        @test isa(constraint_by_name(m, "test2"), Nothing)
        # prepare constraint with duplicate name
        m.constrs[2] = ScalarConstraint(inf + par, MOI.LessThan(1.0))
        m.constr_to_name[2] = "new"
        m.name_to_constr = nothing
        # test for duplciate name error
        @test_throws ErrorException constraint_by_name(m, "new")
    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))
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, x, parameter_bounds = (par == 1))
    data = DiscreteMeasureData(par, [1], [1])
    meas = measure(inf + par - x, data)
    rv = ReducedInfiniteVariableRef(m, 42)
    # test build_constraint (single)
    @testset "JuMP.build_constraint (Single)" begin
        # test bounded constraint
        con = BoundedScalarConstraint(inf, MOI.EqualTo(42.0),
                                ParameterBounds(Dict(par => IntervalSet(0, 1))),
                                ParameterBounds(Dict(par => IntervalSet(0, 1))))
        @test build_constraint(error, inf, MOI.EqualTo(42.0),
             parameter_bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))).func == con.func
        @test build_constraint(error, inf, MOI.EqualTo(42.0),
               parameter_bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))).set == con.set
        @test build_constraint(error, inf, MOI.EqualTo(42.0),
         parameter_bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))).bounds == con.bounds
        # test scalar constraint
        con = ScalarConstraint(inf, MOI.EqualTo(42.0))
        @test build_constraint(error, inf, MOI.EqualTo(42.0)).func == con.func
        @test build_constraint(error, inf, MOI.EqualTo(42.0)).set == con.set
    end
    # test build_constraint (expr)
    @testset "JuMP.build_constraint (Expr)" begin
        # test bounded constraint
        con = BoundedScalarConstraint(inf + pt, MOI.EqualTo(42.0),
                                ParameterBounds(Dict(par => IntervalSet(0, 1))),
                                ParameterBounds(Dict(par => IntervalSet(0, 1))))
        @test build_constraint(error, inf + pt, MOI.EqualTo(42.0),
             parameter_bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))).func == con.func
        @test build_constraint(error, inf + pt, MOI.EqualTo(42.0),
               parameter_bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))).set == con.set
        @test build_constraint(error, inf + pt, MOI.EqualTo(42.0),
         parameter_bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))).bounds == con.bounds
        # test scalar constraint
        con = ScalarConstraint(inf + pt, MOI.EqualTo(42.0))
        @test build_constraint(error, inf + pt, MOI.EqualTo(42.0)).func == con.func
        @test build_constraint(error, inf + pt, MOI.EqualTo(42.0)).set == con.set
    end
    # test _update_var_constr_mapping
    @testset "_update_var_constr_mapping" begin
        # test initial use of variable
        @test isa(InfiniteOpt._update_var_constr_mapping([inf, par, meas, rv], 1),
                  Nothing)
        @test m.var_to_constrs[JuMP.index(inf)] == [1]
        @test m.param_to_constrs[JuMP.index(par)] == [1]
        @test m.meas_to_constrs[JuMP.index(meas)] == [1]
        @test m.reduced_to_constrs[JuMP.index(rv)] == [1]
        # test secondary use of variable
        @test isa(InfiniteOpt._update_var_constr_mapping([inf, par, meas, rv], 2),
                  Nothing)
        @test m.var_to_constrs[JuMP.index(inf)] == [1, 2]
        @test m.param_to_constrs[JuMP.index(par)] == [1, 2]
        @test m.meas_to_constrs[JuMP.index(meas)] == [1, 2]
        @test m.reduced_to_constrs[JuMP.index(rv)] == [1, 2]
        # Undo changes
        delete!(m.var_to_constrs, JuMP.index(inf))
        delete!(m.param_to_constrs, JuMP.index(par))
        delete!(m.meas_to_constrs, JuMP.index(meas))
        delete!(m.reduced_to_constrs, JuMP.index(rv))
    end
    # test _check_and_update_bounds (BoundedScalarConstraint with no hold variables)
    @testset "_check_and_update_bounds (Bounded no Hold)" begin
        c = BoundedScalarConstraint(inf, MOI.EqualTo(0.0), ParameterBounds(),
                                    ParameterBounds())
        @test InfiniteOpt._check_and_update_bounds(m, c, [inf]).bounds == c.bounds
    end
    # test _check_and_update_bounds (BoundedScalarConstraint with hold variables)
    @testset "_check_and_update_bounds (Bounded with Hold)" 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.intervals[par] == IntervalSet(1, 1)
        @test InfiniteOpt._check_and_update_bounds(m, c,
                                      [x, inf]).orig_bounds == ParameterBounds()
    end
    # test _check_and_update_bounds (ScalarConstraint with no hold variables)
    @testset "_check_and_update_bounds (Scalar no Hold)" begin
        c = JuMP.ScalarConstraint(inf, MOI.EqualTo(0.0))
        @test InfiniteOpt._check_and_update_bounds(m, c, [inf]) == c
    end
    # test _check_and_update_bounds (ScalarConstraint with hold variables)
    @testset "_check_and_update_bounds (Scalar with Hold)" begin
        c = JuMP.ScalarConstraint(inf + x, MOI.EqualTo(0.0))
        @test InfiniteOpt._check_and_update_bounds(m, c,
                            [x, inf]).bounds.intervals[par] == IntervalSet(1, 1)
        set_parameter_bounds(x, ParameterBounds(), force = true)
        m.has_hold_bounds = false
    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 parameter only input
        con = ScalarConstraint(par, MOI.EqualTo(42.0))
        @test_throws ErrorException add_constraint(m, con, "a")
        # test rejects bad variables
        @variable(Model(), z)
        con = ScalarConstraint(z, MOI.EqualTo(42.0))
        @test_throws VariableNotOwned add_constraint(m, con, "a")
        @hold_variable(InfiniteModel(), g)
        con = ScalarConstraint(g, MOI.EqualTo(42.0))
        @test_throws VariableNotOwned add_constraint(m, con, "a")
        # test bounded constraint
        index = m.next_constr_index + 1
        con = BoundedScalarConstraint(inf + pt, MOI.EqualTo(42.0),
                                ParameterBounds(Dict(par => IntervalSet(0, 1))),
                                ParameterBounds(Dict(par => IntervalSet(0, 1))))
        @test add_constraint(m, con, "a") == InfiniteConstraintRef(m, index,
                                                                  ScalarShape())
        # test infinite constraint
        con = ScalarConstraint(inf + pt, MOI.EqualTo(42.0))
        @test add_constraint(m, con, "b") == InfiniteConstraintRef(m, index + 1,
                                                                  ScalarShape())
        # test measure constraint
        con = ScalarConstraint(x + meas, MOI.EqualTo(42.0))
        @test add_constraint(m, con, "c") == MeasureConstraintRef(m, index + 2,
                                                                  ScalarShape())
        # test finite constraint
        con = ScalarConstraint(x + pt, MOI.EqualTo(42.0))
        @test add_constraint(m, con, "d") == FiniteConstraintRef(m, index + 3,
                                                                  ScalarShape())
        @test name(FiniteConstraintRef(m, index + 3, ScalarShape())) == "d"
        @test !m.constr_in_var_info[index + 3]
        @test !optimizer_model_ready(m)
        # test with bounded hold variables
        @set_parameter_bounds(x, par == 1)
        con = ScalarConstraint(inf + pt + x, MOI.EqualTo(42.0))
        @test add_constraint(m, con, "b") == InfiniteConstraintRef(m, index + 4,
                                                                  ScalarShape())
        @test m.constrs[index + 4].bounds.intervals[par] == IntervalSet(1, 1)
        @test m.constrs[index + 4].orig_bounds == ParameterBounds()
        set_parameter_bounds(x, ParameterBounds(), force = true)
    end
    # test macro
    @testset "JuMP.@constraint" begin
        # test scalar constraint
        index = m.next_constr_index + 1
        @test @constraint(m, e, x + pt -2 <= 2) == FiniteConstraintRef(m, index,
                                                                  ScalarShape())
        @test isa(m.constrs[index], ScalarConstraint)
        # test bounded scalar constraint
        @test @constraint(m, f, inf + meas <= 2,
            parameter_bounds = ParameterBounds(Dict(par => IntervalSet(0, 1)))) == InfiniteConstraintRef(m, index + 1, ScalarShape())
        @test isa(m.constrs[index + 1], BoundedScalarConstraint)
    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 = :(InfiniteOpt._make_vector(t))
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(Dict($(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
        index = m.next_constr_index + 1
        cref = InfiniteConstraintRef(m, index, ScalarShape())
        @test @BDconstraint(m, par in [0, 1], inf + x == 2) == cref
        @test m.constrs[1].bounds == bounds1
        # test anonymous constraint with comparison
        cref = InfiniteConstraintRef(m, index + 1, ScalarShape())
        @test @BDconstraint(m, 0 <= par <= 1, inf + x == 2) == cref
        @test m.constrs[2].bounds == bounds1
        # test anonymous constraint with equality
        cref = InfiniteConstraintRef(m, index + 2, ScalarShape())
        @test @BDconstraint(m, par == 0, inf + x == 2) == cref
        @test m.constrs[3].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test anonymous multiple constraint
        cref1 = InfiniteConstraintRef(m, index + 3, ScalarShape())
        cref2 = InfiniteConstraintRef(m, index + 4, ScalarShape())
        @test @BDconstraint(m, [1:2](par == 0), inf + x == 2) == [cref1, cref2]
        @test m.constrs[4].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test m.constrs[5].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test anonymous multiple bounds
        cref = InfiniteConstraintRef(m, index + 5, ScalarShape())
        @test @BDconstraint(m, (par == 0, pars[1] in [0, 1]), inf + x == 2) == cref
        @test m.constrs[6].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0),
                                                  pars[1] => IntervalSet(0, 1)))
        # test anonymous multiple bounds with vector input
        cref = InfiniteConstraintRef(m, index + 6, ScalarShape())
        @test @BDconstraint(m, (par == 0, pars in [0, 1]), inf + x == 2) == cref
        @test m.constrs[7].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0),
                    pars[1] => IntervalSet(0, 1), pars[2] => IntervalSet(0, 1)))
        # test named constraint
        cref = InfiniteConstraintRef(m, index + 7, ScalarShape())
        @test @BDconstraint(m, c1(par in [0, 1]), inf + x == 2) == cref
        @test m.constrs[8].bounds == bounds1
        # test named constraint with other vector type
        cref = InfiniteConstraintRef(m, index + 8, ScalarShape())
        @test @BDconstraint(m, c2(par in [0, 1], 0 <= pars <= 1), inf + x == 2) == cref
        @test m.constrs[9].bounds == ParameterBounds(Dict(par => IntervalSet(0, 1),
                    pars[1] => IntervalSet(0, 1), pars[2] => IntervalSet(0, 1)))
        # test named constraints
        cref1 = InfiniteConstraintRef(m, index + 9, ScalarShape())
        cref2 = InfiniteConstraintRef(m, index + 10, ScalarShape())
        @test @BDconstraint(m, c3[1:2](par == 0), inf + x == 2) == [cref1, cref2]
        @test m.constrs[10].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test m.constrs[11].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test different container type
        cref1 = InfiniteConstraintRef(m, index + 11, ScalarShape())
        cref2 = InfiniteConstraintRef(m, index + 12, ScalarShape())
        expected = convert(JuMPC.SparseAxisArray, [cref1, cref2])
        @test @BDconstraint(m, c4[1:2](par == 0), inf + x == 2,
                            container = SparseAxisArray) == expected
        @test m.constrs[12].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test m.constrs[13].bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        # test with SparseAxisArray input
        cref = InfiniteConstraintRef(m, index + 13, ScalarShape())
        pars = convert(JuMPC.SparseAxisArray, pars)
        @test @BDconstraint(m, c5(pars in [0, 1]), inf + x == 2) == cref
        @test m.constrs[14].bounds == ParameterBounds(Dict(pars[1] => IntervalSet(0, 1),
                                                  pars[2] => IntervalSet(0, 1)))
        # test with hold variable bounds
        @set_parameter_bounds(x, pars[1] in [0, 2])
        cref = InfiniteConstraintRef(m, index + 14, ScalarShape())
        @test @BDconstraint(m, c10(par in [0, 1]), inf + x == 2) == cref
        @test m.constrs[15].bounds == ParameterBounds(Dict(par => IntervalSet(0, 1),
                                                  pars[1] => IntervalSet(0, 2)))
        @test m.constrs[15].orig_bounds == ParameterBounds(Dict(par => IntervalSet(0, 1)))
    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 has_parameter_bounds
    @testset "has_parameter_bounds" begin
        @test has_parameter_bounds(c1)
        @test !has_parameter_bounds(c2)
    end
    # test parameter_bounds
    @testset "parameter_bounds" begin
        @test_throws ErrorException parameter_bounds(c2)
        @test parameter_bounds(c1) == ParameterBounds(Dict(par => IntervalSet(0, 1)))
    end
    # 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_bound" begin
        # test already has bounds
        @test isa(add_parameter_bound(c1, par, 0, 1), Nothing)
        @test parameter_bounds(c1) == ParameterBounds(Dict(par => IntervalSet(0, 1)))
        # test doesn't have bounds
        @constraint(m, c3, inf + pt == 0)
        @test isa(add_parameter_bound(c3, par, 0, 1), Nothing)
        @test parameter_bounds(c3) == ParameterBounds(Dict(par => IntervalSet(0, 1)))
    end
    # test delete_parameter_bound
    @testset "delete_parameter_bound" begin
        # test partial deletion
        @test isa(add_parameter_bound(c1, pars[1], 0, 1), Nothing)
        @add_parameter_bounds(x, par == 0)
        @test isa(delete_parameter_bound(c1, par), Nothing)
        expected = ParameterBounds(Dict(par => IntervalSet(0, 0),
                                   pars[1] => IntervalSet(0, 1)))
        @test parameter_bounds(c1) == expected
        # test again without hold bounds
        delete_parameter_bounds(x)
        expected = ParameterBounds(Dict(pars[1] => IntervalSet(0, 1)))
        @test parameter_bounds(c1) == expected
        # test already gone
        @test isa(delete_parameter_bound(c1, par), Nothing)
        @test parameter_bounds(c1) == expected
    end
    # test delete_parameter_bounds
    @testset "delete_parameter_bounds" begin
        # test partial deletion
        @test isa(add_parameter_bound(c1, pars[1], 0, 1), 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
    # 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[2] == 0)
        @test isa(@set_parameter_bounds(c1, par == 1, force = true), Nothing)
        @test parameter_bounds(c1) == ParameterBounds(Dict(par => IntervalSet(1, 1),
                                                  pars[2] => 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).intervals[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).intervals[pars[2]] == IntervalSet(0, 0)
        @test parameter_bounds(c2).intervals[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, HoldVariableRef, MOI.LessThan) == 1
        @test num_constraints(m, InfiniteVariableRef, MOI.GreaterThan) == 1
    end
    # test search function type
    @testset "Function" begin
        @test num_constraints(m, HoldVariableRef) == 3
        @test num_constraints(m, InfiniteVariableRef) == 1
    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 = [FiniteConstraintRef(m, 4, ScalarShape())]
        @test all_constraints(m, HoldVariableRef, MOI.LessThan) == list
        list = [InfiniteConstraintRef(m, 1, ScalarShape())]
        @test all_constraints(m, InfiniteVariableRef, MOI.GreaterThan) == list
    end
    # test search function type
    @testset "Function" begin
        list = [FiniteConstraintRef(m, 3, ScalarShape()),
                FiniteConstraintRef(m, 4, ScalarShape()),
                FiniteConstraintRef(m, 5, ScalarShape())]
        @test all_constraints(m, HoldVariableRef) == list
        list = [InfiniteConstraintRef(m, 1, ScalarShape())]
        @test all_constraints(m, InfiniteVariableRef) == list
    end
    # test search set type
    @testset "Set" begin
        list = [FiniteConstraintRef(m, 4, ScalarShape())]
        @test all_constraints(m, MOI.LessThan) == list
        list = [InfiniteConstraintRef(m, 1, ScalarShape()),
                FiniteConstraintRef(m, 2, ScalarShape()),
                FiniteConstraintRef(m, 3, ScalarShape())]
        @test all_constraints(m, MOI.GreaterThan) == list
    end
    # test search total
    @testset "Total" begin
        list = [InfiniteConstraintRef(m, 1, ScalarShape()),
                FiniteConstraintRef(m, 2, ScalarShape()),
                FiniteConstraintRef(m, 3, ScalarShape()),
                FiniteConstraintRef(m, 4, ScalarShape()),
                FiniteConstraintRef(m, 5, ScalarShape())]
        @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 = [(InfiniteVariableRef, MOI.GreaterThan{Float64}),
                (PointVariableRef, MOI.GreaterThan{Float64}),
                (HoldVariableRef, MOI.GreaterThan{Float64}),
                (HoldVariableRef, MOI.LessThan{Float64}),
                (HoldVariableRef, MOI.Integer)]
    @test list_of_constraint_types(m) == expected
end


text:
module HorseML

using LinearAlgebra
using Statistics

export Regression, Classification, Preprocessing, LossFunction, Tree, NeuralNetwork, Clustering, gpu, cpu, fit!

include("utils.jl")
include("Regression/Regression.jl")
include("Classification/Classification.jl")
include("Tree/Tree.jl")
include("NeuralNetwork/NN.jl")
include("Clustering/Clustering.jl")
include("Preprocessing/Preprocessing.jl")
include("LossFunction/LossFunctions.jl")
include("gpu.jl")

end

ext: jl


text:
function mark_cross!(r::Robot)
    a = 0
    b = 0
    while (isborder(r,Nord) && isborder(r,West))==false
        while !isborder(r , Nord)
            move!(r, Nord)
            a = a + 1
        end
        while !isborder(r, West)
            move!(r,West)
            b = b + 1
        end
    end
    putmarker!(r)
    while !isborder(r, Ost)
        move!(r, Ost)
    end
    putmarker!(r)
    while !isborder(r, Sud)
        move!(r, Sud)
    end
    putmarker!(r)
    while !isborder(r, West)
        move!(r, West)
    end
    putmarker!(r)
    while !isborder(r, Nord)
        move!(r, Nord)
    end
    while b > 0
        move!(r, Ost)
        b = b - 1
    end
    while a > 0
        if isborder(r,Sud)
            s = 0
            while isborder(r, Sud)
                move!(r, West)
                s = s + 1
            end
            move!(r, Sud)
            a = a - 1
            while isborder(r, Ost)
                move!(r, Sud)
                a = a - 1
            end
            while s > 0
                move!(r, Ost)
                s = s - 1
            end
        else
            move!(r, Sud)
            a = a - 1
    end
end


text:
module githubactions4julia

# Write your package code here.
"""
githubactions4juliasample()
return true
"""
function githubactions4juliasample()
    return true
end

export githubactions4juliasample
end


text:
using Zygote

begin
W = rand(2, 3); x = rand(3);
@show gradient(W -> sum(W*x), W)
# Since we only have W as parameter, the first element is the grad
# a 2x3 array
@show gradient(W -> sum(W*x), W)[1]
end

begin
d = Dict()
# Note the do is a x-> ....  function
# Gradient of that function is calculated at x=5
# function acts like x**2  so  its gradient is 2x=2*5=10
@show  gradient(5) do x
         d[:x] = x
         d[:x] * d[:x]
       end

@show  d[:x]
end

max_stars_repo_name: daryoush/CudaAndFluxExperiments

max_issues_repo_name: daryoush/CudaAndFluxExperiments

max_forks_repo_name: daryoush/CudaAndFluxExperiments


text:
@testset "Infinite ODE" begin
    """
    ODEsolver reports `NaN` (None) when steady state is absent.
    """

    function infinite!(du,u,p,t)
        du[1] = (u[1]-0.5)^2 * p[1] + 0.5 # Impossible to be zero
    end

    de = DEsteady(
        func = infinite!,
        u0 = [0.4],
        p = [1.0],
        method=SSRootfind()
    )

    sol = solve(de)

    @test all(isnan.(sol.u)) #NaN
end


text:
using Interpolations
tsolparallel1

range(1.0, stop = 400.0, length = 2000)

lambda_k = 1:0.2:1000

lambda_omega = range(1.0, stop = 2000.0, length = 4000)

f(x1, x2) = log(x1 + x2)

function getImlambdaData(k, omega)
    n = round(Int64, (omega - 1) / step(lambda_omega)) + 1
    tsolparallel1[n](k)[1]
end
function getRelambdaData(k, omega)
    n = round(Int64, (omega - 1) / step(lambda_omega)) + 1
    tsolparallel1[n](k)[2]
end

ImA = [getImlambdaData(x1, x2) for x1 in lambda_k, x2 in lambda_omega]
ReA = [getRelambdaData(x1, x2) for x1 in lambda_k, x2 in lambda_omega]

Imitp = interpolate(ImA, BSpline(Linear()))
Imsitp = scale(Imitp, lambda_k, lambda_omega)
Imetpf = extrapolate(Imsitp, Flat())

Reitp = interpolate(ReA, BSpline(Linear()))
Resitp = scale(Reitp, lambda_k, lambda_omega)
Reetpf = extrapolate(Resitp, Flat())

plot(k -> -getImlambdaData(k, lambda_omega[800]), 1, 800)

plot!(k -> -Imetpf(k, lambda_omega[1000] + 100), 1, 800)

m2fun(x) = sol1(x)[1]

twopointflow(100.0, 10.0, m2fun(10.0), 145.0, Imetpf)

using HCubature

function Imgamma2(p0)
    -hcubature(
        x -> twopointflow(p0, x[1], m2fun(x[1]), 145.0, Imetpf),
        [1.0],
        [800.0],
        atol = 1e-4,
        rtol = 1e-4,
    )[1]
end

function Regamma2(p0)
    -hcubature(
        x -> twopointflow(p0, x[1], m2fun(x[1]), 145.0, Reetpf),
        [1.0],
        [800.0],
        atol = 1e-4,
        rtol = 1e-4,
    )[1] - p0^2 + m2fun(800.0)
end

function bosonicspec(p0)
    -2 * Imgamma2(p0) / (Imgamma2(p0)^2 + Regamma2(p0)^2)
end

plot(p0 -> Imgamma2(p0), 1.0, 400.0)
plot(p0 -> Regamma2(p0), 1.0, 400.0)

plot(p0 -> -bosonicspec(p0), 1.0, 400.0, yaxis = :log)


text:
struct NullData end
const DEFAULT_DATA = Iterators.cycle((NullData(),))
Base.iterate(::NullData, i=1) = nothing
Base.length(::NullData) = 0

get_maxiters(data) = Iterators.IteratorSize(typeof(DEFAULT_DATA)) isa Iterators.IsInfinite ||
                     Iterators.IteratorSize(typeof(DEFAULT_DATA)) isa Iterators.SizeUnknown ?
                     typemax(Int) : length(data)

function update!(x::AbstractArray, x̄::AbstractArray{<:ForwardDiff.Dual})
  x .-= x̄
end

function update!(x::AbstractArray, x̄)
  x .-= getindex.(ForwardDiff.partials.(x̄),1)
end

function update!(opt, x, x̄)
  x .-= Flux.Optimise.apply!(opt, x, x̄)
end

function update!(opt, x, x̄::AbstractArray{<:ForwardDiff.Dual})
  x .-= Flux.Optimise.apply!(opt, x, getindex.(ForwardDiff.partials.(x̄),1))
end

function update!(opt, xs::Flux.Zygote.Params, gs)
  for x in xs
    gs[x] == nothing && continue
    update!(opt, x, gs[x])
  end
end

abstract type DiffMode end
struct ForwardDiffMode <: DiffMode end
struct TrackerDiffMode <: DiffMode end
struct ReverseDiffMode <: DiffMode end
struct ZygoteDiffMode <: DiffMode end

macro withprogress(progress, exprs...)
  quote
    if $progress
      $DiffEqBase.maybe_with_logger($DiffEqBase.default_logger($Logging.current_logger())) do
        $ProgressLogging.@withprogress $(exprs...)
      end
    else
      $(exprs[end])
    end
  end |> esc
end

"""
    sciml_train(loss, θ, opt, data = DEFAULT_DATA; cb, maxiters)

Optimizes the `loss(θ,curdata...)` function with respect to the parameter vector
`θ` iterating over the `data`. By default the data iterator is empty, i.e.
`loss(θ)` is used. The first output of the loss function is considered the loss.
Extra outputs are passed to the callback.

The keyword arguments are as follows:

- `cb`: For specifying a callback function `cb(θ,l...)` which acts on the current
  parameters `θ`, where `l...` are the returns of the latest call to `loss`. This
  callback should return a boolean where, if true, the training will end prematurely.
- `maxiters`: Specifies the maximum number of iterations for the optimization.
  Required if a Flux optimizer is chosen and no data iterator is given, otherwise
  defaults to infinite.
- `save_best`: Specifies whether you want the best solution, corresponding to the
  lowest loss function value. If false, the last solution found will be returned.
"""
function sciml_train(loss, _θ, opt, _data = DEFAULT_DATA;
                     cb = (args...) -> false,
                     maxiters = get_maxiters(data),
                     progress=true, save_best=true)

  # Flux is silly and doesn't have an abstract type on its optimizers, so assume
  # this is a Flux optimizer
  θ = copy(_θ)
  ps = Flux.params(θ)

  if _data == DEFAULT_DATA && maxiters == typemax(Int)
    error("For Flux optimizers, either a data iterator must be provided or the `maxiters` keyword argument must be set.")
  elseif _data == DEFAULT_DATA && maxiters != typemax(Int)
    data = Iterators.repeated((), maxiters)
  elseif maxiters != typemax(Int)
    data = take(_data, maxiters)
  else
    data = _data
  end

  t0 = time()

  local x, min_err
  min_err = typemax(eltype(_θ)) #dummy variables
  min_opt = 1

  @withprogress progress name="Training" begin
    for (i,d) in enumerate(data)
      gs = Flux.Zygote.gradient(ps) do
        x = loss(θ,d...)
        first(x)
      end
      cb_call = cb(θ,x...)
      if !(typeof(cb_call) <: Bool)
        error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
      elseif cb_call
        break
      end
      msg = @sprintf("loss: %.3g", x[1])
      progress && ProgressLogging.@logprogress msg i/maxiters
      update!(opt, ps, gs)

      if save_best
        if first(x) < first(min_err)  #found a better solution
          min_opt = opt
          min_err = x
        end
        if i == maxiters  #Last iteration, revert to best.
          opt = min_opt
          cb(θ,min_err...)
        end
      end
    end
  end

  _time = time()

  Optim.MultivariateOptimizationResults(opt,
                                        _θ,# initial_x,
                                        θ, #pick_best_x(f_incr_pick, state),
                                        first(x), # pick_best_f(f_incr_pick, state, d),
                                        maxiters, #iteration,
                                        maxiters >= maxiters, #iteration == options.iterations,
                                        true, # x_converged,
                                        0.0,#T(options.x_tol),
                                        0.0,#T(options.x_tol),
                                        NaN,# x_abschange(state),
                                        NaN,# x_abschange(state),
                                        true,# f_converged,
                                        0.0,#T(options.f_tol),
                                        0.0,#T(options.f_tol),
                                        NaN,#f_abschange(d, state),
                                        NaN,#f_abschange(d, state),
                                        true,#g_converged,
                                        0.0,#T(options.g_tol),
                                        NaN,#g_residual(d),
                                        false, #f_increased,
                                        nothing,
                                        maxiters,
                                        maxiters,
                                        0,
                                        true,
                                        NaN,
                                        _time-t0)
end

decompose_trace(trace::Optim.OptimizationTrace) = last(trace)
decompose_trace(opt::BlackBoxOptim.OptRunController) = BlackBoxOptim.best_candidate(opt)
decompose_trace(trace) = trace

function sciml_train(loss, θ, opt::Optim.AbstractOptimizer, data = DEFAULT_DATA;
                      cb = (args...) -> (false), maxiters = get_maxiters(data),
                      diffmode = ZygoteDiffMode(), kwargs...)
  local x, cur, state
  cur,state = iterate(data)

  function _cb(trace)
    cb_call = cb(decompose_trace(trace).metadata["x"],x...)
    if !(typeof(cb_call) <: Bool)
      error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
    end
    cur,state = iterate(data,state)
    cb_call
  end

  _loss = function (θ)
    x = loss(θ,cur...)
    first(x)
  end

  if diffmode isa ZygoteDiffMode
    f!  =  function (θ)
      _loss(θ)
    end

    g!  =  function (G,θ)
      _x,lambda = Flux.Zygote.pullback(_loss,θ)
      grad = first(lambda(1))
      if eltype(grad) <: ForwardDiff.Dual
        G .= getindex.(ForwardDiff.partials.(grad),1)
      else
        G .= grad
      end
      return _x
    end

    fg! =  function (G,θ)
      _x,lambda = Flux.Zygote.pullback(_loss,θ)
      if G != nothing
        grad = first(lambda(1))
        if eltype(grad) <: ForwardDiff.Dual
          G .= getindex.(ForwardDiff.partials.(grad),1)
        else
          G .= grad
        end
      end
      return _x
    end

    if opt isa Optim.KrylovTrustRegion
      hv! = function (H,θ,v)
        _θ = ForwardDiff.Dual.(θ,v)
        H .= getindex.(ForwardDiff.partials.(Flux.Zygote.gradient(_loss,_θ)[1]),1)
      end
      optim_f = Optim.TwiceDifferentiableHV(f!,fg!,hv!,θ)
    else
      h! = function (H,θ)
        H .= ForwardDiff.jacobian(θ) do θ
          Flux.Zygote.gradient(_loss,θ)[1]
        end
      end
      optim_f = TwiceDifferentiable(f!, g!, fg!, h!, θ)
    end

  else
    if diffmode isa TrackerDiffMode
      optim_fg! =  function (F,G,θ)
        _x,lambda = Tracker.forward(_loss,θ)
        if G != nothing
          G .= first(lambda(1))
        end
        if F != nothing
          return _x
        end
        return _x
      end
    elseif diffmode isa ForwardDiffMode
      optim_fg! = let res = DiffResults.GradientResult(θ)
        function (F,G,θ)
          if G != nothing
            ForwardDiff.gradient!(res,_loss,θ)
            G .= DiffResults.gradient(res)
          end

          if F != nothing
            return DiffResults.value(res)
          end

          return _x
        end
      end
    elseif diffmode isa ReverseDiffMode
      optim_fg! = let res = DiffResults.GradientResult(θ)
        function (F,G,θ)
          if G != nothing
            ReverseDiff.gradient!(res,_loss,θ)
            G .= DiffResults.gradient(res)
          end

          if F != nothing
            return DiffResults.value(res)
          end

          return _x
        end
      end
    end
    optim_f = Optim.only_fg!(optim_fg!)
  end
  Optim.optimize(optim_f, θ, opt,
           Optim.Options(;extended_trace=true,callback = _cb,
                         f_calls_limit = maxiters, kwargs...))
end

function sciml_train(loss, θ, opt::Optim.AbstractConstrainedOptimizer,
                     data = DEFAULT_DATA;
                     lower_bounds, upper_bounds,
                     cb = (args...) -> (false), maxiters = get_maxiters(data))
  local x, cur, state
  cur,state = iterate(data)

  function _cb(trace)
    cb_call = cb(decompose_trace(trace).metadata["x"],x...)
    if !(typeof(cb_call) <: Bool)
      error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
    end
    cur,state = iterate(data,state)
    cb_call
  end

  function optim_fg!(F,G,θ)
    _x,lambda = Flux.Zygote.pullback(θ) do θ
      x = loss(θ,cur...)
      first(x)
    end

    if G != nothing
      G .= first(lambda(1))
    end

    if F != nothing
      return _x
    end

    return _x
  end

  Optim.optimize(Optim.only_fg!(optim_fg!), lower_bounds, upper_bounds, θ, opt,
           Optim.Options(extended_trace=true,callback = _cb,
                         f_calls_limit = maxiters))
end

struct BBO
  method::Symbol
end

BBO() = BBO(:adaptive_de_rand_1_bin)

function sciml_train(loss, _θ, opt::BBO = BBO(), data = DEFAULT_DATA;lower_bounds, upper_bounds,
                      cb = BlackBoxOptim.trace_progress, maxiters = get_maxiters(data), kwargs...)
  local x, cur, state
  cur,state = iterate(data)

  function _cb(trace)
    cb_call = cb(decompose_trace(trace),x...)
    if !(typeof(cb_call) <: Bool)
      error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
    end
    if cb_call == true
      BlackBoxOptim.shutdown_optimizer!(trace) #doesn't work
    end
    cur,state = iterate(data,state)
    cb_call
  end

  _loss = function (θ)
    x = loss(θ,cur...)
    first(x)
  end

  bboptre = BlackBoxOptim.bboptimize(_loss;Method = opt.method, SearchRange = [(lower_bounds[i], upper_bounds[i]) for i in 1:length(lower_bounds)], MaxSteps = maxiters, CallbackFunction = _cb, CallbackInterval = 0.0, kwargs...)

  Optim.MultivariateOptimizationResults(opt.method,
                                        [NaN],# initial_x,
                                        best_candidate(bboptre), #pick_best_x(f_incr_pick, state),
                                        best_fitness(bboptre), # pick_best_f(f_incr_pick, state, d),
                                        bboptre.iterations, #iteration,
                                        bboptre.iterations >= maxiters, #iteration == options.iterations,
                                        true, # x_converged,
                                        0.0,#T(options.x_tol),
                                        0.0,#T(options.x_tol),
                                        NaN,# x_abschange(state),
                                        NaN,# x_abschange(state),
                                        true,# f_converged,
                                        0.0,#T(options.f_tol),
                                        0.0,#T(options.f_tol),
                                        NaN,#f_abschange(d, state),
                                        NaN,#f_abschange(d, state),
                                        true,#g_converged,
                                        0.0,#T(options.g_tol),
                                        NaN,#g_residual(d),
                                        false, #f_increased,
                                        nothing,
                                        maxiters,
                                        maxiters,
                                        0,
                                        true,
                                        NaN,
                                        bboptre.elapsed_time)
end


text:
function _compute_sum_rows_cols(A)
    R_mx1 = zeros(eltype(A), size(A, 1))
    C_nx1 = zeros(eltype(A), size(A, 2))

    sum_row = A * ones( size(A, 2) )
    EPS = eps(eltype(A)) # choose float32 to make it more sparse
    R_mx1[sum_row .> EPS] .= 1 ./ (sum_row[sum_row .> EPS])

    sum_col = A' * ones( size(A, 1) )
    C_nx1[sum_col .> EPS] .= 1 ./ (sum_col[sum_col .> EPS])

    return R_mx1, C_nx1
end

@doc raw"""
    recon2d_sirt

Reconstruct a 2d image by SIRT

# Args
A : Forward opeartor
b : Projection data
u0: Initial guess of image
niter: number of iterations

u <- u + CA'R(b - Au)
"""
function recon2d_sirt!(u0::Array{T, 2}, A::SparseMatrixCSC{T,Int}, b0::Array{T, 2}, niter::Int; min_value=nothing, max_value=nothing) where {T <: AbstractFloat}
    R_mx1, C_nx1 = _compute_sum_rows_cols(A)

    At = sparse(A') # this significatnly improves the performance
    u = vec(u0)
    b = vec(b0)
    r = similar(b)

    for it = 1:niter
        # println("$it")
        r .= b .- A*u
        u .+= C_nx1 .* (At * (R_mx1 .* r))

        if !isnothing(min_value) # much faster than min_value == nothing
            u .= max.(u, min_value)
        end
        if !isnothing(max_value)
            u .= min.(u, max_value)
        end
        if it % 30 == 0
            residual = sum( r .^ 2 )  / length(r)
            println("$it l2 residual: $residual")
        end
    end

    residual = sum( r .^ 2 )  / length(r)
end

# ---------------------------------------------------------
# ---------------------------------------------------------
# ---------------------------------------------------------

@doc raw"""
    recon2d_slices_sirt

Reconstruct 3d image by SIRT slice by slice

# Args
A : Forward opeartor
b : Projection data [H x W x nangles]
u0: Initial guess of 3D image [H x W x nslice]
niter: number of iterations

u <- u + CA'R(b - Au)
"""
function recon2d_stack_sirt!(u::Array{T, 3}, A::SparseMatrixCSC{T,Int}, b_::Array{T, 3}, niter::Int; min_value=nothing) where {T <: AbstractFloat}
    if size(b_, 1) != size(u, 3)
        error("Not supported yet when the detector row count == image slice size")
    end

    R_mx1, C_nx1 = _compute_sum_rows_cols(A)
    At = A' # this significatnly improves the performance
    # RA = R_mx1 .* A
    # CAt = C_nx1 .* A'

    nslice = size(u, 3)

    b_axWxH = permutedims(b_, [3, 2, 1])
    bR = reshape(b_axWxH, :, nslice)
    println(size(bR))

    halfslice = Int(floor(nslice/2))
    H, W = size(u, 1), size(u, 2)

    u_view = reshape(u, :, nslice)
    u_temp = similar(u_view)
    temp = zeros(size(bR))

    use_multi_threads = true
    println("@ preprocessing done, run SIRT.")
    println("@ num of julia threads: $(Threads.nthreads())")
    # multi threaded version
    for it = 1:niter
        Threads.@threads for slice=1:nslice
            bR_ = view(bR, :, slice)
            uu = view(u_view, :, slice)
            uu_temp = view(u_temp, :, slice)
            temp_slice = view(temp, :, slice)

            temp_slice .= R_mx1 .* (bR_ .- mul!(temp_slice, A, uu))
            uu .+= C_nx1 .* mul!(uu_temp, At, temp_slice)
        end

        if !isnothing(min_value)
            u .= max.(u, min_value)
        end

        if it % 10 == 0
            # residual = sum( temp .^ 2 )  / length(temp)
            println("iter: $it, scaled residual $(sum(temp .^ 2) / length(temp))")
        end
    end

    residual = (sum(temp .^ 2) / length(temp))
end

# temp .= R_mx1 .* (b .- mul!(temp, A, u_view ) )
# mul!(u_view, At, temp)
# u_view .*= C_nx1

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
#------------------------------------------------------------------------------

function common_denominator(coeffs::Vector{Rational{BigInt}})
    den = BigInt(1)
    for c in coeffs
        # TODO: Modular arithmetic
        den = lcm(den, denominator(c))
    end
    den
end

# given an array of rational coeffs
# scales them by the minimal number such that the coeffs are all integer
function remove_content!(coeffs::Vector{Rational{BigInt}})
    intcoeffs = Vector{BigInt}(undef, length(coeffs))
    # lcm of map(denominator, coeffs[i])
    # so that den × coeffs[i] is integer
    den = common_denominator(coeffs)
    for i in 1:length(coeffs)
        # TODO: MutableArithmetics
        delta = divexact(den, denominator(coeffs[i]))
        intcoeffs[i] = numerator(coeffs[i]) * delta
    end
    intcoeffs
end

function scale_denominators!(coeffs::Vector{Rational{BigInt}})
    remove_content!(coeffs)
end

# scales denominators inplace, returning new vector of integer coefficients
function scale_denominators!(coeffs::Vector{Vector{Rational{BigInt}}})
    intcoeffs = Vector{Vector{BigInt}}(undef, length(coeffs))
    for i in 1:length(coeffs)
        intcoeffs[i] = scale_denominators!(coeffs[i])
    end
    intcoeffs
end

#------------------------------------------------------------------------------

# coerce each coeff in `coeffs` into residuals modulo `prime`
function reduce_modulo(coeffs::Vector{BigInt}, prime::Int)
    ffcoeffs = Vector{UInt64}(undef, length(coeffs))
    for i in 1:length(coeffs)
        # TODO: faster
        # convert is guaranteed to be exact
        c = coeffs[i]

        if c < 0
            c = c - prime*(div(c, prime) - 1)
        end
        ffcoeffs[i] = UInt64(c % prime)
    end
    ffcoeffs
end

# coerce each coeff in `coeffs` into residuals modulo `prime`
function reduce_modulo(
        coeffs::Vector{Vector{BigInt}},
        ring::PolyRing, prime::Int)

    ffcoeffs = Vector{Vector{UInt64}}(undef, length(coeffs))
    for i in 1:length(coeffs)
        ffcoeffs[i] = reduce_modulo(coeffs[i], prime)
    end
    ring_ff = PolyRing(ring.nvars, ring.explen, ring.ord, prime, :abstract)
    ring_ff, ffcoeffs
end

#------------------------------------------------------------------------------

# Not used
# checks if coefficients in coeffs are not large compared to modulo
# Having relatively large coeffs would mean that reconstruction probably failed
function reconstruction_check(
        coeffs::Vector{Vector{Rational{BigInt}}},
        modulo::BigInt
    )

    @inbounds for i in 1:length(coeffs)
        for j in 1:length(coeffs[i])
            c = coeffs[i][j]
            # heuristic
            # TODO: a better one
            pval = ( numerator(c) + denominator(c) ) ^ 2
            if pval >= modulo
                return false
            end
        end
    end

    return true
end

#------------------------------------------------------------------------------

function isluckyprime(
        coeffs::Vector{Vector{BigInt}},
        prime::Int)
    for poly in coeffs
        for c in poly
            if c % prime == 0
                 return false
            end
        end
    end
    return true
end

function nextluckyprime(
        coeffs::Vector{Vector{BigInt}},
        prevprime::Int)

    if prevprime == 1
        candidate = 2^31 - 1
    else
        candidate = nextprime(prevprime + 1)
    end
    @assert candidate isa Int

    if !isluckyprime(coeffs, candidate)
        candidate = nextluckyprime(coeffs, candidate)
    end

    candidate
end

#------------------------------------------------------------------------------

function nextgoodprime(
        coeffs::Vector{Vector{BigInt}},
        moduli::Vector{Int},
        prevprime::Int)

    p = nextluckyprime(coeffs, prevprime)
    if p in moduli
        return nextgoodprime(coeffs, moduli, p)
    end
    return p
end