text:
"""
    test_quadratic_integration(model::MOI.ModelLike, config::Config)

Test the problem:
```
Min x^2 + xy + y^2 + yz + z^2
st  x + 2y + 3z >= 4 (c1)
    x +  y      >= 1 (c2)
    x, y, z in R
```
"""
function test_quadratic_integration(model::MOI.ModelLike, config::Config)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Float64},
        MOI.GreaterThan{Float64},
    )
    v = MOI.add_variables(model, 3)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3
    cf1 =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], v), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.GreaterThan(4.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ) == 1
    end
    c2 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [v[1], v[2]]),
            0.0,
        ),
        MOI.GreaterThan(1.0),
    )
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ) == 2
    end
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    obj = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.(
            [2.0, 1.0, 2.0, 1.0, 2.0],
            v[[1, 1, 2, 2, 3]],
            v[[1, 2, 2, 3, 3]],
        ),
        MOI.ScalarAffineTerm{Float64}[],
        0.0,
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        obj,
    )
    if _supports(config, MOI.ObjectiveFunction)
        @test obj ≈ MOI.get(
            model,
            MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        )
        @test cf1 ≈ MOI.get(model, MOI.ConstraintFunction(), c1)
        @test MOI.GreaterThan(4.0) == MOI.get(model, MOI.ConstraintSet(), c1)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 13 / 7, config)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), v),
            [4 / 7, 3 / 7, 6 / 7],
            config,
        )
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            # The dual constraint gives
            # [λ_c1 + λ_c2, 2λ_c1 + λ_c2, 3λ_c1] = [2 1 0; 1 2 1; 0 1 2] * x
            #                                    = [11, 16, 15] / 7
            # hence λ_c1 = 5/7 and λ_c2 = 6/7.
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c1),
                5 / 7,
                config,
            )
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c2),
                6 / 7,
                config,
            )
        end
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_integration),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [5 / 7, 6 / 7],
        ),
    )
    return
end

"""
    test_quadratic_duplicate_terms(model::MOI.ModelLike, config::Config)

Test the problem:
```
Min x^2 + xy + y^2 + yz + z^2
st  x + 2y + 3z >= 4 (c1)
    x +  y      >= 1 (c2)
    x, y, z in R
```
"""
function test_quadratic_duplicate_terms(model::MOI.ModelLike, config::Config)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Float64},
        MOI.GreaterThan{Float64},
    )
    v = MOI.add_variables(model, 3)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3
    c1f =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], v), 0.0)
    c1 = MOI.add_constraint(model, c1f, MOI.GreaterThan(4.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ) == 1
    end
    c2 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [v[1], v[2]]),
            0.0,
        ),
        MOI.GreaterThan(1.0),
    )
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ) == 2
    end
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    obj = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.(
            [2.0, 0.5, 0.5, 2.0, 1.0, 1.0, 1.0],
            [v[1], v[1], v[1], v[2], v[2], v[3], v[3]],
            [v[1], v[2], v[2], v[2], v[3], v[3], v[3]],
        ),
        MOI.ScalarAffineTerm.(0.0, v),
        0.0,
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        obj,
    )
    if _supports(config, MOI.ObjectiveFunction)
        @test obj ≈ MOI.get(
            model,
            MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        )
        @test c1f ≈ MOI.get(model, MOI.ConstraintFunction(), c1)
        @test MOI.GreaterThan(4.0) == MOI.get(model, MOI.ConstraintSet(), c1)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 13 / 7, config)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), v),
            [4 / 7, 3 / 7, 6 / 7],
            config,
        )
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c1),
                5 / 7,
                config,
            )
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c2),
                6 / 7,
                config,
            )
        end
    end
    # change objective to Max -2(x^2 + xy + y^2 + yz + z^2)
    # First clear the objective, this is needed if a `Bridges.Objective.SlackBridge` is used.
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    obj2 = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.(
            [-4.0, -1.0, -1.0, -4.0, -2.0, -2.0, -2.0],
            [v[1], v[1], v[1], v[2], v[2], v[3], v[3]],
            [v[1], v[2], v[2], v[2], v[3], v[3], v[3]],
        ),
        MOI.ScalarAffineTerm.(0.0, v),
        0.0,
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        obj2,
    )
    if _supports(config, MOI.ObjectiveFunction)
        @test obj2 ≈ MOI.get(
            model,
            MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        )
    end
    if _supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(
            MOI.get(model, MOI.ObjectiveValue()),
            -2 * 13 / 7,
            config,
        )
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), v),
            [4 / 7, 3 / 7, 6 / 7],
            config,
        )
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c1),
                10 / 7,
                config,
            )
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c2),
                12 / 7,
                config,
            )
        end
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_duplicate_terms),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [5 / 7, 6 / 7],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [10 / 7, 12 / 7],
        ),
    )
    return
end

"""
    test_quadratic_nonhomogeneous(model::MOI.ModelLike, config::Config)

Test the problem:
```
minimize 2 x^2 + y^2 + xy + x + y + 1
    s.t.  x, y >= 0
          x + y = 1
```
"""
function test_quadratic_nonhomogeneous(model::MOI.ModelLike, config::Config)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.VariableIndex,
        MOI.GreaterThan{Float64},
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Float64},
        MOI.EqualTo{Float64},
    )
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    c1 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
            0.0,
        ),
        MOI.EqualTo(1.0),
    )
    vc1 = MOI.add_constraint(model, x, MOI.GreaterThan(0.0))
    # We test this after the creation of every `VariableIndex` constraint
    # to ensure a good coverage of corner cases.
    @test vc1.value == x.value
    vc2 = MOI.add_constraint(model, y, MOI.GreaterThan(0.0))
    @test vc2.value == y.value
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    obj = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([4.0, 2.0, 1.0], [x, y, x], [x, y, y]),
        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
        1.0,
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        obj,
    )
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 2.875, config)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), [x, y]),
            [0.25, 0.75],
            config,
        )
        @test isapprox(
            MOI.get(model, MOI.ConstraintPrimal(), vc1),
            0.25,
            config,
        )
        @test isapprox(
            MOI.get(model, MOI.ConstraintPrimal(), vc2),
            0.75,
            config,
        )
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            # The dual constraint gives
            # [λ_vc1 + λ_c1, λ_vc2 + λ_c1] = [4 1; 1 2] * x + [1, 1] = [11, 11] / 4
            # since `vc1` and `vc2` are not active, `λ_vc1` and `λ_vc2` are
            # zero so `λ_c1 = 11/4`.
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c1),
                11 / 4,
                config,
            )
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), vc1), 0, config)
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), vc2), 0, config)
        end
    end
    # change back to linear
    #        max 2x + y + 1
    #       s.t.  x, y >= 0
    #             x + y = 1
    # (x,y) = (1,0), obj = 3
    # First clear the objective, this is needed if a
    # `Bridges.Objective.SlackBridge` is used.
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    objf =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], [x, y]), 1.0)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        objf,
    )
    if _supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 3.0, config)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), [x, y]),
            [1.0, 0.0],
            config,
        )
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), vc1), 1.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), vc2), 0.0, config)
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), c1), -2, config)
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), vc1), 0, config)
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), vc2), 1, config)
        end
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_nonhomogeneous),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1 / 4, 3 / 4],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [11 / 4],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 0.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [-2.0],
        ),
    )
    return
end

"""
    test_quadratic_constraint_integration(model::MOI.ModelLike, config::Config)

Test the problem:
```
Max x  + y
st -x  + y >= 0 (c1[1])
    x  + y >= 0 (c1[2])
    x² + y <= 2 (c2)
Optimal solution
x = 1/2, y = 7/4
```
"""
function test_quadratic_constraint_integration(
    model::MOI.ModelLike,
    config::Config,
)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Float64},
        MOI.Nonnegatives,
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        MOI.LessThan{Float64},
    )
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    c1 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 1, 2, 2],
                MOI.ScalarAffineTerm.([-1.0, 1.0, 1.0, 1.0], [x, y, x, y]),
            ),
            [0.0, 0.0],
        ),
        MOI.Nonnegatives(2),
    )
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.VectorAffineFunction{Float64},
                MOI.Nonnegatives,
            }(),
        ) == 1
    end
    c2f = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(2.0, x, x)],
        [MOI.ScalarAffineTerm(1.0, y)],
        0.0,
    )
    c2 = MOI.add_constraint(model, c2f, MOI.LessThan(2.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ) == 1
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    if _supports(config, MOI.ConstraintFunction)
        @test isapprox(
            MOI.get(model, MOI.ConstraintFunction(), c2),
            c2f,
            config,
        )
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 2.25, config)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), [x, y]),
            [0.5, 1.75],
            config,
        )
        @test isapprox(
            MOI.get(model, MOI.ConstraintPrimal(), c1),
            [5 / 4, 9 / 4],
            config,
        )
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), c2), 2, config)
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            # The dual constraint gives
            # [-1, -1] = [-λ_c1[1] + λ_c1[2], λ_c1[1] + λ_c1[2]] + λ_c2 * ([2 0; 0 0] * x + [0, 1])
            #          = [-λ_c1[1] + λ_c1[2] + λ_c2, λ_c1[1] + λ_c1[2] + λ_c2]
            # By complementary slackness, we have `λ_c1 = [0, 0]` so
            # `λ_c2 = -1`.
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c1),
                [0, 0],
                config,
            )
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), c2), -1, config)
        end
    end
    # try delete quadratic constraint and go back to linear
    # FIXME why is this commented ?
    # MOI.delete(model, c2)
    #
    # MOI.optimize!(model)
    #
    # @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
    #
    # @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    #
    # @test MOI.get(model, MOI.ObjectiveValue()), 0.0 atol=atol rtol=rtol
    return
end

function setup_test(
    ::typeof(test_quadratic_constraint_integration),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1 / 2, 7 / 4],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [zeros(2)],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
    )
    return
end

"""
    test_quadratic_constraint_basic(model::MOI.ModelLike, config::Config)

Test the problem:
```
 Max x
s.t. x^2 <= 2 (c)
```
"""
function test_quadratic_constraint_basic(model::MOI.ModelLike, config::Config)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        MOI.LessThan{Float64},
    )
    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    cf = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(2.0, x, x)],
        [MOI.ScalarAffineTerm(0.0, x)],
        0.0,
    )
    c = MOI.add_constraint(model, cf, MOI.LessThan(2.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ) == 1
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    if _supports(config, MOI.ConstraintFunction)
        @test isapprox(MOI.get(model, MOI.ConstraintFunction(), c), cf, config)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), √2, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), x), √2, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), c), 2, config)
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            # The dual constraint gives
            # -1 = λ * (2 * x) = λ * 2√2
            # hence λ = -1/(2√2).
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c),
                -1 / (2 * √2),
                config,
            )
        end
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_constraint_basic),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [√2],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1 / (2 * √2)],
        ),
    )
    return
end

"""
    test_quadratic_constraint_minimize(
        model::MOI.ModelLike,
        config::Config,
    )

Test the problem:
```
 Min -x
s.t.  x^2 <= 2
```
"""
function test_quadratic_constraint_minimize(
    model::MOI.ModelLike,
    config::Config,
)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        MOI.LessThan{Float64},
    )
    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    cf = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(2.0, x, x)],
        MOI.ScalarAffineTerm{Float64}[],
        0.0,
    )
    c = MOI.add_constraint(model, cf, MOI.LessThan(2.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ) == 1
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    if _supports(config, MOI.ConstraintFunction)
        @test isapprox(MOI.get(model, MOI.ConstraintFunction(), c), cf, config)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
        end
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), -√2, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), x), √2, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), c), 2, config)
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            # The dual constraint gives
            # -1 = λ * (2 * x) = λ * 2√2
            # hence λ = -1/(2√2).
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c),
                -1 / (2 * √2),
                config,
            )
        end
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_constraint_minimize),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [√2],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1 / (2 * √2)],
        ),
    )
    return
end

"""
Max  x
s.t. x^2 + x * y + y^2 <= 3 (c)
     y == 1
"""
function _test_quadratic_constraint_helper(
    model::MOI.ModelLike,
    config::Config,
    less_than::Bool,
)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    quad_set = less_than ? MOI.LessThan(3.0) : MOI.GreaterThan(-3.0)
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        typeof(quad_set),
    )
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    vc = MOI.add_constraint(model, y, MOI.EqualTo(1.0))
    @test vc.value == y.value
    cf = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([2.0, 1.0, 2.0], [x, x, y], [x, y, y]),
        [MOI.ScalarAffineTerm(0.0, x)],
        0.0,
    )
    if !less_than
        MOIU.operate!(-, Float64, cf)
    end
    c = MOI.add_constraint(model, cf, quad_set)
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarQuadraticFunction{Float64},
                typeof(quad_set),
            }(),
        ) == 1
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    if _supports(config, MOI.ConstraintFunction)
        @test isapprox(MOI.get(model, MOI.ConstraintFunction(), c), cf, config)
        @test quad_set == MOI.get(model, MOI.ConstraintSet(), c)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 1.0, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), x), 1.0, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), y), 1.0, config)
        @test isapprox(
            MOI.get(model, MOI.ConstraintPrimal(), c),
            MOI.constant(quad_set),
            config,
        )
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), vc), 1.0, config)
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            # The dual constraint gives
            # [-1, 0] = λ_c * [2 1; 1 2] * x + [0, λ_vc]
            #         = [3λ_c, 3λ_c + λ_vc]
            # hence `λ_c = -1/3` and `λ_vc = 1`.
            λ_c = less_than ? -1 / 3 : 1 / 3
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), c), λ_c, config)
            @test isapprox(MOI.get(model, MOI.ConstraintDual(), vc), 1, config)
        end
    end
    return
end

function test_quadratic_constraint_LessThan(
    model::MOI.ModelLike,
    config::Config,
)
    _test_quadratic_constraint_helper(model, config, true)
    return
end

function setup_test(
    ::typeof(test_quadratic_constraint_LessThan),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 1.0],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1 / 3],
        ),
    )
    return
end

function test_quadratic_constraint_GreaterThan(
    model::MOI.ModelLike,
    config::Config,
)
    _test_quadratic_constraint_helper(model, config, false)
    return
end

function setup_test(
    ::typeof(test_quadratic_constraint_GreaterThan),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 1.0],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}) => [1 / 3],
        ),
    )
    return
end

"""
    test_quadratic_nonconvex_constraint_integration(
        model::MOI.ModelLike,
        config::Config,
    )

Test the problem:
```
Max 2x + y
s.t. x * y <= 4 (c)
     x, y >= 1
```
"""
function test_quadratic_nonconvex_constraint_integration(
    model::MOI.ModelLike,
    config::Config,
)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        MOI.LessThan{Float64},
    )
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    vc1 = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    @test vc1.value == x.value
    vc2 = MOI.add_constraint(model, y, MOI.GreaterThan(1.0))
    @test vc2.value == y.value
    cf = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(1.0, x, y)],
        [MOI.ScalarAffineTerm(0.0, x)],
        0.0,
    )
    c = MOI.add_constraint(model, cf, MOI.LessThan(4.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ) == 1
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([2.0, 1.0], [x, y]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    if _supports(config, MOI.ConstraintFunction)
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 9.0, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), x), 4.0, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), y), 1.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), c), 4.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), vc1), 4.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), vc2), 1.0, config)
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_nonconvex_constraint_integration),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, [4.0, 1.0], MOI.FEASIBLE_POINT),
    )
    return
end

"""
    test_quadratic_nonconvex_constraint_basic(
        model::MOI.ModelLike,
        config::Config,
    )

Test the problem:
```
Find x, y
s.t. x * y == 4 (c)
     x * x == 4 (c2)
     x, y >= 0
```
"""
function test_quadratic_nonconvex_constraint_basic(
    model::MOI.ModelLike,
    config::Config,
)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        MOI.EqualTo{Float64},
    )
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    vc1 = MOI.add_constraint(model, x, MOI.GreaterThan(0.0))
    @test vc1.value == x.value
    vc2 = MOI.add_constraint(model, y, MOI.GreaterThan(0.0))
    @test vc2.value == y.value
    cf = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(1.0, x, y)],
        [MOI.ScalarAffineTerm(0.0, x)],
        0.0,
    )
    c = MOI.add_constraint(model, cf, MOI.EqualTo(4.0))
    cf2 = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(2.0, x, x)],
        [MOI.ScalarAffineTerm(0.0, x)],
        0.0,
    )
    c2 = MOI.add_constraint(model, cf2, MOI.EqualTo(4.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.EqualTo{Float64},
            }(),
        ) == 2
    end
    if _supports(config, MOI.ConstraintFunction)
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
        @test cf2 ≈ MOI.get(model, MOI.ConstraintFunction(), c2)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), x), 2.0, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), y), 2.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), c), 4.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), c2), 4.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), vc1), 2.0, config)
        @test isapprox(MOI.get(model, MOI.ConstraintPrimal(), vc2), 2.0, config)
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_nonconvex_constraint_basic),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, [2.0, 2.0], MOI.FEASIBLE_POINT),
    )
    return
end

"""
    test_quadratic_SecondOrderCone_basic(model::MOI.ModelLike, config::Config)

Test the non-convex formulation of a SecondOrderCone constraint:
```
min t
s.t. x + y >= 1 (c1)
     x^2 + y^2 <= t^2 (c2)
     t >= 0 (bound)
```

This is a useful test because some solvers have support for detecting that it is
actually a convex problem due to the `t >= 0` bound.

If your solver does not have this support (and therefore fails this test), you
can safely exclude it.
"""
function test_quadratic_SecondOrderCone_basic(
    model::MOI.ModelLike,
    config::Config,
)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        MOI.LessThan{Float64},
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Float64},
        MOI.GreaterThan{Float64},
    )
    @requires MOI.supports_constraint(
        model,
        MOI.VariableIndex,
        MOI.GreaterThan{Float64},
    )
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    t = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3
    c1f =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
    c1 = MOI.add_constraint(model, c1f, MOI.GreaterThan(1.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ) == 1
    end
    c2f = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([2.0, 2.0, -2.0], [x, y, t], [x, y, t]),
        MOI.ScalarAffineTerm{Float64}[],
        0.0,
    )
    c2 = MOI.add_constraint(model, c2f, MOI.LessThan(0.0))
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ) == 1
    end
    bound = MOI.add_constraint(model, t, MOI.GreaterThan(0.0))
    @test bound.value == t.value
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.VariableIndex,
                MOI.GreaterThan{Float64},
            }(),
        ) == 1
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, t)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    if _supports(config, MOI.ConstraintFunction)
        @test c1f ≈ MOI.get(model, MOI.ConstraintFunction(), c1)
        @test c2f ≈ MOI.get(model, MOI.ConstraintFunction(), c2)
    end
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 1 / √2, config)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), [x, y, t]),
            [0.5, 0.5, 1 / √2],
            config,
        )
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), [t, x, y, t]),
            [1 / √2, 0.5, 0.5, 1 / √2],
            config,
        )
        if _supports(config, MOI.ConstraintDual)
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
            # The dual constraint gives
            # [1, 0, 0] = [λ_bound, λ_c1, λ_c1] + λ_c2 * [-2 0 0; 0 2 0; 0 0 2] * x
            #           = [λ_bound - √2*λ_c2, λ_c1 + λ_c2, λ_c1 + λ_c2]
            # and since `bound` is not active, `λ_bound = 0`.
            # This gives `λ_c2 = -1/√2` and `λ_c1 = 1/√2`.
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c1),
                1 / √2,
                config,
            )
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), c2),
                -1 / √2,
                config,
            )
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), bound),
                0,
                config,
            )
        end
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_SecondOrderCone_basic),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1 / 2, 1 / 2, 1 / √2],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1 / √2],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1 / √2],
        ),
    )
    return
end

"""
    test_quadratic_Integer_SecondOrderCone(model::MOI.ModelLike, config::Config)

Run an integration test on the problem:
```
 min    - 2y - 1z
s.t.   x           == 1
      [x,  y,   z] in SecondOrderCone()
           y,   z  in ZeroOne()
```
"""
function test_quadratic_Integer_SecondOrderCone(
    model::MOI.ModelLike,
    config::Config,
)
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Float64},
        MOI.Zeros,
    )
    @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @requires MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.SecondOrderCone,
    )
    x, y, z = MOI.add_variables(model, 3)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-2.0, -1.0], [y, z]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
            [-1.0],
        ),
        MOI.Zeros(1),
    )
    MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x, y, z]),
        MOI.SecondOrderCone(3),
    )
    if _supports(config, MOI.NumberOfConstraints)
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ) == 1
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SecondOrderCone}(),
        ) == 1
    end
    loc = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test length(loc) == 2
    @test (MOI.VectorAffineFunction{Float64}, MOI.Zeros) in loc
    @test (MOI.VectorOfVariables, MOI.SecondOrderCone) in loc
    bin1 = MOI.add_constraint(model, y, MOI.ZeroOne())
    # We test this after the creation of every `VariableIndex` constraint
    # to ensure a good coverage of corner cases.
    @test bin1.value == y.value
    bin2 = MOI.add_constraint(model, z, MOI.ZeroOne())
    @test bin2.value == z.value
    if _supports(config, MOI.optimize!)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), -2, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), x), 1, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), y), 1, config)
        @test isapprox(MOI.get(model, MOI.VariablePrimal(), z), 0, config)
    end
    return
end

function setup_test(
    ::typeof(test_quadratic_Integer_SecondOrderCone),
    model::MOI.Utilities.MockOptimizer,
    config::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, [1.0, 1.0, 0.0]),
    )
    return
end


text:
export findfirst, findall
################################################################
# Find Node
################################################################
# findfirst
function Base.findfirst(name::String, doc::Document, argAmlType::Type{<:AbsNode})
    elm = findfirst(name, root(doc), argAmlType)
    return elm
end

function Base.findfirst(name::String, node::Node, argAmlType::Type{<:AbsNormal})
    # if hasdocument(node)
    #     elm = findfirst(name, node)
    # else
        elm = findfirstelm(name, node)
    # end
    return elm
end

function Base.findfirst(name::String, node::Node, argAmlType::Type{AbsAttribute})
    elms = findfirstatt(name, node)
    return elms
end

function Base.findfirst(indexstr::String, node::Node, argAmlType::Type{AbsText})
    index = parse_textindex(indexstr)
    elm = findfirsttext(index, node)
    return elm
end
################################################################
# findall
function Base.findall(name::String, doc::Document, argAmlType::Type{<:AbsNode})
    elms = findall(name, root(doc), argAmlType)
    return elms
end

function Base.findall(name::String, node::Node, argAmlType::Type{<:AbsNormal})
    # if hasdocument(node)
    #     elms = findall(name, node) # a vector of Node elements
    # else
        elms = findallelm(name, node) # a vector of Node elements
    # end
    return elms
end

function Base.findall(name::String, node::Node, argAmlType::Type{AbsAttribute})
    elms = findallatt(name, node)
    return elms
end

function Base.findall(indicesstr::String, node::Node, argAmlType::Type{AbsText})
    indices = parse_textindices(indicesstr)
    elms = findalltext(indices, node)
    return elms
end

################################################################
# Utilities
################################################################
function Base.get(node::Node, name::String, defval)
    if haskey(node, name)
        return node[name]
    else
        return defval
    end
end
################################################################

"""
    findfirstatt(name, node)

Finds the first attribute with `name`

findfirst with ignoring namespaces. It considers att.name for returning the elements

Much faster than EzXML.findfirst
"""
function findfirstatt(name::String, node::Node)
    out = nothing # return nothing if nothing is found
    for att in eachattribute(node)
        if att.name == name
            out = att
            break
        end
    end
    return out
end

"""
    findallatt(string, node)

Finds all the attributes with `name`

findallatt with ignoring namespaces. It considers att.name for returning the elements

Much faster than EzXML.findall
"""
function findallatt(name::String, node::Node)
    out = Node[]
    for att in eachattribute(node)
        if att.name == name
            push!(out, att)
        end
    end
    if !isempty(out)
        return out
    else # return nothing if nothing is found
        return nothing
    end
end

################################################################
# Local searchers (no namespace)
"""
    findfirstelm(string, node)

findfirst with ignoring namespaces. It considers element.name for returning the elements

Much faster than EzXML.findfirst
"""
function findfirstelm(name::String, node::Node)
    out = nothing # return nothing if nothing is found
    for child in eachelement(node)
        if child.name == name
            out = child
            break
        end
    end
    return out
end

"""
    findallelm(string, node)

findallelm with ignoring namespaces. It considers element.name for returning the elements

Much faster than EzXML.findall
"""
function findallelm(name::String, node::Node)
    out = Node[]
    for child in eachelement(node)
        if child.name == name
            push!(out, child)
        end
    end
    if !isempty(out)
        return out
    else # return nothing if nothing is found
        return nothing
    end
end

################################################################
# Text node utils

function findtext(indexstr::String, node::Node)
    if indexstr == ""
        index = 1
    else
        index = eval(Meta.parse(indexstr))
    end
    xpath = "text()[position() = $index]"
    out = findfirst(xpath, node)
    return out
end

# function findvecttext(indexstr::String, node::Node)
#     if indexstr == ""
#         index = Colon()
#     else
#         index = eval(Meta.parse(indexstr))
#     end
#     xpath = "text()"
#     out = findfirst(xpath, node)
#     return out
# end

"""
    parse_textindex(indexstr::String)

Index is a String of an Integer e.g. "2". If indexstr is empty (`""`) it returns the first one found.
"""
function parse_textindex(indexstr::String)
    if indexstr == ""
        index = 1
    elseif indexstr == "end"
        index = Inf
    else
        indexExpr = Meta.parse(indexstr)
        indexExpr isa Integer || error("give index as an Integer e.g. \"2\"")
        index = eval(indexExpr)
    end
    return index
end

"""
    findfirsttext(index::Integer, node)

finds the text node at position given by index.

faster than `findtext()`
"""
function findfirsttext(index::Integer, node::Node)
    iText = 0
    out = nothing # return nothing if nothing is found
    for child in eachnode(node)
        if istext(child)
            iText +=1
            if iText == index
                out = child
                break
            end
        end
    end
    return out
end
function findfirsttext(index::Float64, node::Node)
    if index != Inf
        error("index should be \"end\"")
    end
    out = nothing # return nothing if nothing is found
    for child in eachnode(node)
        if istext(child)
            out = child
        end
    end
    return out
end

"""
    parse_textindices(indicesstr::String)

Index is a String of an Integer e.g. "2". If indicesstr is empty (`""`) it returns the all of the text nodes.
"""
function parse_textindices(indicesstr::String)
    if indicesstr == ""
        indices = Colon()
    else
        indicesExpr = Meta.parse(indicesstr)
        indicesExpr.head == :vect || error("give indices as a vetor e.g. [2:3], [2, 3] ,[:]")
        indices = eval(indicesExpr)
        # TODO find a better way
        if typeof(indices) <: Vector{<:Union{UnitRange{Int64},Colon}}
            indices = indices[1]
        end
    end
    return indices
end

"""
    findalltext(indices, node)

finds the text node at positions given by indices.

faster than `findvecttext()`
"""
function findalltext(indices::Colon, node::Node)
    out = Node[]
    for child in eachnode(node)
        if istext(child)
            push!(out, child)
        end
    end
    if !isempty(out)
        return out
    else # return nothing if nothing is found
        return nothing
    end
end

function findalltext(indices::AbstractVector, node::Node)
    out = Vector{Node}(undef, length(indices))
    iTextLocation = 1
    outIndex = 1
    for child in eachnode(node)
        if istext(child)
            if iTextLocation in indices
                out[outIndex] = child
            end
            outIndex += 1
        else
            iTextLocation +=1
        end
    end
    if !isempty(out)
        return out
    else # return nothing if nothing is found
        return nothing
    end
end

################################################################

ext: jl


text:
using Nabla

# Generate the values required to compute a matrix quadratic form.
N = 5
B = randn(5, 5)
A = B.'B + UniformScaling(1e-6)

# Low-level API computation of derivatives.

# Construct some trackable data.
x_, y_ = randn(N), randn(N)
x, y = Leaf.(Tape(), (x_, y_))

# Compute the forward pass.
z = x.' * (A * y) # Temporary bracketting because we don't support RowVectors yet.

println("Output of the forward pass is:")
println(z)
println()
println("y is $(Nabla.unbox(z)).")
println()

# Get the reverse tape.
z̄ = ∇(z)
println("Output of reverse-pass is")
println(z̄)
println()

# Index into the reverse tape using x to get the gradient of `y` w.r.t. `x`.
x̄ = z̄[x]
println("Gradient of z w.r.t. x at $x_ is $x̄.")
println()

ȳ = z̄[y]
println("Gradient of z w.r.t. y at $y_ is $ȳ")

# (Current) High-Level API computation of derivatives. I will probably maintain this
# interface and extend is significantly as it is currently rather limited. It just returns
# a function which returns the gradient, which isn't really what you want.

# Define the function to be differentiated. Parameters w.r.t. which we want gradients must
# be arguments. Parameters that we don't want gradients w.r.t. should be passed in via a
# closure.
@unionise f(x::AbstractVector, y::AbstractVector) = x.'A * y

# Compute a function `∇f` which computes the derivative of `f` w.r.t. the inputs.
∇f = ∇(f)

# Compute the derivative of `f` w.r.t. `x` at `x_`. Result is currently a 1-Tuple. Might
# introduce a special case for unary functions where it just returns the result.
(x̄, ȳ) = ∇f(x_, y_)

@assert x̄ == z̄[x]
@assert ȳ == z̄[y]


text:
"""
Cache for all model results
"""
mutable struct ResultCache
    data::Dict{ParamCacheKey, ParamResultCache}
    max_size::Int
    min_flush_size::Int
end

function ResultCache()
    return ResultCache(Dict{ParamCacheKey, ParamResultCache}(), 0, 0)
end

# PERF: incremental improvement if we manually keep track of the current size.
# Could be prone to bugs if we miss a change.
# The number of containers is not expected to be more than 100.

#decrement_size!(cache::ResultCache, size) = cache.size -= size
#increment_size!(cache::ResultCache, size) = cache.size += size

get_max_size!(cache::ResultCache) = cache.max_size
get_min_flush_size(cache::ResultCache) = cache.min_flush_size
get_param_cache(cache::ResultCache, key) = cache.data[key]
get_size(cache::ResultCache) = reduce(+, (get_size(x) for x in values(cache.data)))
is_full(cache::ResultCache, cur_size) = cur_size >= cache.max_size
set_max_size!(cache::ResultCache, x) = cache.max_size = x
set_min_flush_size!(cache::ResultCache, x) = cache.min_flush_size = x

function add_param_cache!(cache::ResultCache, key, flush_rule)
    cache.data[key] = ParamResultCache(key, flush_rule)
    #@debug "Added cache container for" key flush_rule
end

"""
Return true if the cache has data that has not been flushed to storage.
"""
function has_dirty(cache::ResultCache)
    for param_cache in values(cache.data)
        if has_dirty(param_cache)
            return true
        end
    end

    return false
end

"""
Return true if the data for `timestamp` is stored in cache.
"""
function is_cached(cache::ResultCache, type, name, problem, timestamp)
    return is_cached(cache, make_cache_key(problem, type, name), timestamp)
end

is_cached(cache::ResultCache, key, timestamp) = has_timestamp(cache.data[key], timestamp)

"""
Log the cache hit percentages for all caches.
"""
function log_cache_hit_percentages(cache::ResultCache)
    for key in sort!(collect(keys(cache.data)))
        param_cache = cache.data[key]
        cache_hit_pecentage = get_cache_hit_percentage(param_cache)
        @debug "Cache stats" key cache_hit_pecentage
    end
end

"""
Read the result from cache. Callers must first call [`is_cached`](@ref) to check if the
timestamp is present.
"""
function read_result(cache::ResultCache, problem, type, name, timestamp)
    return read_result(cache, make_cache_key(problem, type, name), timestamp)
end

read_result(cache::ResultCache, key, timestamp) = cache.data[key].data[timestamp]


text:
# include("./oneBody/group.jl")
# using .OneBody

function markov(state, groups, steps, updates, weight = ones(length(updates)), rng = Random.GLOBAL_RNG)
    z = sum(weight)
    N = length(updates)
    prop = [sum(weight[1:iw]) / z for (iw, w) in enumerate(weight)]
    println(prop)
    for s in 1:steps
        x = rand(rng)
        for (i, p) in prop
            if x < p
                # if updates(i)
                # end
            end
        end
        return state, groups
    end

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
@doc raw"""
    kf = KalmanFilter(x, F, H, B, P, Q, R, offset)

Create a KalmanFilter object.

for the dynamical system
```math
\begin{aligned}
x_{t+1} & = Fx_t + Bu_t + w_t + \text{offset}\\
y_t     & = Hx_t + v_t,
\end{aligned}
```
where $x_t\in \mathbb{R}^n$ are the states, $u\in \mathbb{R}^p$ are known inputs
 $y_t\in \mathbb{R}^m$ are measured outputs and
 $w_t\in \mathbb{R}^n,\ v_t \in \mathbb{R}^m$ are unknown, uncorrelated disturbances,
 with positive definite covariance matrices $Q \in \mathbb{R}^{n\times n}$ and
 $R \in \mathbb{R}^{m\times m}$ respectively. The offset is an $n$-dimensional constant.
"""
function KalmanFilter(
    x::AbstractVector{T},
    F::AbstractMatrix{T},
    H::AbstractMatrix{T},
    B::AbstractMatrix{T},
    P::AbstractMatrix{T},
    Q::AbstractMatrix{T},
    R::AbstractMatrix{T},
    offset::AbstractVector{T},
    ) where T<: AbstractFloat

    n = size(x)[1]
    m = size(H)[1]
    p = size(B)[2]
    # Validate Args
    size(F)[1] != n ? DimensionMismatch("Dimension mismatch of F") |> throw : nothing
    size(H)[2] != n ? DimensionMismatch("Dimension mismatch of H") |> throw : nothing
    size(B)[1] != n ? DimensionMismatch("Dimension mismatch of B") |> throw : nothing
    size(P)[1] != n ? DimensionMismatch("P and F incompatible") |> throw : nothing
    size(Q)[1] != n ? DimensionMismatch("Q and F incompatible") |> throw : nothing
    size(R)[1] != m ? DimensionMismatch("R and C incompatible") |> throw : nothing
    !isposdef(P) ? DomainError("P must be positive-definite") |> throw : nothing
    !isposdef(Q) ? DomainError("Q must be positive-definite") |> throw : nothing
    !isposdef(R) ? DomainError("R must be positive-definite") |> throw : nothing

    c = Base.RefValue{T}(0)
    pdy = Base.RefValue{T}(1)
    x0 = copy(x)
    F = copy(F)
    H = copy(H)
    B = copy(B)
    P = copy(P)
    Q = copy(Q)
    R = copy(R)
    offset = copy(offset)
    KalmanFilter{T}(x0, F, H, B, P, Q, R, offset, c, pdy)
end

# lazy wrappers
function KalmanFilter(x0, F, H, P0, Q, R)
    size(x0) == () ? x0 = [x0] : nothing
    B = zeros(eltype(x0),size(x0,1))
    KalmanFilter(x0, F, H, B, P0, Q, R)
end

function KalmanFilter(x0, F, H, B, P0, Q, R)
    size(x0) == () ? x0 = [x0] : nothing
    offset = zeros(eltype(x0),size(x0))
    KalmanFilter(x0, F, H, B, P0, Q, R, offset)
end

function KalmanFilter(x0, F, H, B, P0, Q, R, offset)
    size(x0) == () ? x0 = [x0] : nothing
    size(F) == () ? F = fill(F,1,1) : nothing
    size(H) == () ? H = fill(H,1,1) : nothing
    size(B) == () ? B = fill(B,1,1) :
        typeof(B) <: AbstractVector ? B = reshape(B, (:,1)) : nothing
    size(P0) == () ? P0 = fill(P0, 1,1) : nothing
    size(Q) == () ? Q = fill(Q, 1,1) : nothing
    size(R) == () ? R = fill(R, 1,1) : nothing
    size(offset) == () ? offset = [offset] : nothing

    x0, F, H, B, P0, Q, R, offset = promote_arrays(x0, F, H, B, P0, Q, R, offset)

    return KalmanFilter(x0, F, H, B, P0, Q, R, offset)

end

# Utility functions
function promote_arrays(arrays...)
    eltype = Base.promote_eltype(arrays...)
    eltype <: AbstractFloat ? nothing : eltype = Float64
    return tuple([convert(Array{eltype}, array) for array in arrays]...)
end


text:
module AbstractDifferentiation

using LinearAlgebra, ExprTools, Requires, Compat

export AD

const AD = AbstractDifferentiation

abstract type AbstractBackend end
abstract type AbstractFiniteDifference <: AbstractBackend end
abstract type AbstractForwardMode <: AbstractBackend end
abstract type AbstractReverseMode <: AbstractBackend end

struct HigherOrderBackend{B} <: AbstractBackend
    backends::B
end
reduce_order(b::AbstractBackend) = b
function reduce_order(b::HigherOrderBackend)
    if length(b.backends)==1
        return lowest(b) # prevent zero tuple and subsequent error when reducing over HigherOrderBackend
    else
        return HigherOrderBackend(reverse(Base.tail(reverse(b.backends))))
    end
end
lowest(b::AbstractBackend) = b
lowest(b::HigherOrderBackend) = b.backends[end]
second_lowest(b::AbstractBackend) = b
second_lowest(b::HigherOrderBackend) = lowest(reduce_order(b))

# If the primal value is in y, extract it.
# Otherwise, re-compute it, e.g. in finite diff.
primal_value(::AbstractFiniteDifference, ::Any, f, xs) = f(xs...)
primal_value(::AbstractBackend, ys, ::Any, ::Any) = primal_value(ys)
primal_value(x::Tuple) = map(primal_value, x)
primal_value(x) = x

function derivative(ab::AbstractBackend, f, xs::Number...)
    der = getindex.(jacobian(lowest(ab), f, xs...), 1)
    if der isa Tuple
        return der
    else
        return (der,)
    end
end

function gradient(ab::AbstractBackend, f, xs...)
    return reshape.(adjoint.(jacobian(lowest(ab), f, xs...)),size.(xs))
end
function jacobian(ab::AbstractBackend, f, xs...) end
function jacobian(ab::HigherOrderBackend, f, xs...)
    jacobian(lowest(ab), f, xs...)
end

function hessian(ab::AbstractBackend, f, x)
    if x isa Tuple
        # only support computation of Hessian for functions with single input argument
        @assert length(x) == 1
        x = x[1]
    end
    return jacobian(second_lowest(ab), x -> begin
        gradient(lowest(ab), f, x)[1] # gradient returns a tuple
    end, x)
end

function value_and_derivative(ab::AbstractBackend, f, xs::Number...)
    value, jacs = value_and_jacobian(lowest(ab), f, xs...)
    return value[1], getindex.(jacs, 1)
end
function value_and_gradient(ab::AbstractBackend, f, xs...)
    value, jacs = value_and_jacobian(lowest(ab), f, xs...)
    return value, reshape.(adjoint.(jacs),size.(xs))
end
function value_and_jacobian(ab::AbstractBackend, f, xs...)
    local value
    primalcalled = false
    if lowest(ab) isa AbstractFiniteDifference
        value = primal_value(ab, nothing, f, xs)
        primalcalled = true
    end
    jacs = jacobian(lowest(ab), (_xs...,) -> begin
        v = f(_xs...)
        if !primalcalled
            value = primal_value(ab, v, f, xs)
            primalcalled = true
        end
        return v
    end, xs...)

    return value, jacs
end
function value_and_hessian(ab::AbstractBackend, f, x)
    if x isa Tuple
        # only support computation of Hessian for functions with single input argument
        @assert length(x) == 1
        x = x[1]
    end

    local value
    primalcalled = false
    if ab isa AbstractFiniteDifference
        value = primal_value(ab, nothing, f, (x,))
        primalcalled = true
    end
    hess = jacobian(second_lowest(ab), _x -> begin
        v, g = value_and_gradient(lowest(ab), f, _x)
        if !primalcalled
            value = primal_value(ab, v, f, (x,))
            primalcalled = true
        end
        return g[1] # gradient returns a tuple
    end, x)
    return value, hess
end
function value_and_hessian(ab::HigherOrderBackend, f, x)
    if x isa Tuple
        # only support computation of Hessian for functions with single input argument
        @assert length(x) == 1
        x = x[1]
    end
    local value
    primalcalled = false
    hess = jacobian(second_lowest(ab), (_x,) -> begin
        v, g = value_and_gradient(lowest(ab), f, _x)
        if !primalcalled
            value = primal_value(ab, v, f, (x,))
            primalcalled = true
        end
        return g[1]  # gradient returns a tuple
    end, x)
    return value, hess
end
function value_gradient_and_hessian(ab::AbstractBackend, f, x)
    if x isa Tuple
        # only support computation of Hessian for functions with single input argument
        @assert length(x) == 1
        x = x[1]
    end
    local value
    primalcalled = false
    grads, hess = value_and_jacobian(second_lowest(ab), _x -> begin
        v, g = value_and_gradient(lowest(ab), f, _x)
        if !primalcalled
            value = primal_value(second_lowest(ab), v, f, (x,))
            primalcalled = true
        end
        return g[1] # gradient returns a tuple
    end, x)
    return value, (grads,), hess
end
function value_gradient_and_hessian(ab::HigherOrderBackend, f, x)
    if x isa Tuple
        # only support computation of Hessian for functions with single input argument
        @assert length(x) == 1
        x = x[1]
    end
    local value
    primalcalled = false
    grads, hess = value_and_jacobian(second_lowest(ab), _x -> begin
        v, g = value_and_gradient(lowest(ab), f, _x)
        if !primalcalled
            value = primal_value(second_lowest(ab), v, f, (x,))
            primalcalled = true
        end
        return g[1] # gradient returns a tuple
    end, x)
    return value, (grads,), hess
end

function pushforward_function(
    ab::AbstractBackend,
    f,
    xs...,
)
    return (ds) -> begin
        return jacobian(lowest(ab), (xds...,) -> begin
            if ds isa Tuple
                @assert length(xs) == length(ds)
                newxs = xs .+ ds .* xds
                return f(newxs...)
            else
                @assert length(xs) == length(xds) == 1
                newx = xs[1] + ds * xds[1]
                return f(newx)
            end
        end, _zero.(xs, ds)...)
    end
end
function value_and_pushforward_function(
    ab::AbstractBackend,
    f,
    xs...,
)
    return (ds) -> begin
        if !(ds isa Tuple)
            ds = (ds,)
        end
        @assert length(ds) == length(xs)
        local value
        primalcalled = false
        if ab isa AbstractFiniteDifference
            value = primal_value(ab, nothing, f, xs)
            primalcalled = true
        end
        pf = pushforward_function(lowest(ab), (_xs...,) -> begin
            vs = f(_xs...)
            if !primalcalled
                value = primal_value(lowest(ab), vs, f, xs)
                primalcalled = true
            end
            return vs
        end, xs...)(ds)

        return value, pf
    end
end

_zero(::Number, d::Number) = zero(d)
_zero(::Number, d::AbstractVector) = zero(d)
_zero(::AbstractVector, d::AbstractVector) = zero(eltype(d))
_zero(::AbstractVector, d::AbstractMatrix) = zero(similar(d, size(d, 2)))
_zero(::AbstractMatrix, d::AbstractMatrix) = zero(d)
_zero(::Any, d::Any) = zero(d)

@inline _dot(x, y) = dot(x, y)
@inline function _dot(x::AbstractVector, y::UniformScaling)
    @assert length(x) == 1
    return @inbounds dot(x[1], y.λ)
end
@inline function _dot(x::AbstractVector, y::AbstractMatrix)
    @assert size(y, 2) == 1
    return dot(x, y)
end

function pullback_function(ab::AbstractBackend, f, xs...)
    return (ws) -> begin
        return gradient(lowest(ab), (xs...,) -> begin
            vs = f(xs...)
            if ws isa Tuple
                @assert length(vs) == length(ws)
                return sum(Base.splat(_dot), zip(ws, vs))
            else
                return _dot(vs, ws)
            end
        end, xs...)
    end
end
function value_and_pullback_function(
    ab::AbstractBackend,
    f,
    xs...,
)
    return (ws) -> begin
        local value
        primalcalled = false
        if ab isa AbstractFiniteDifference
            value = primal_value(ab, nothing, f, xs)
            primalcalled = true
        end
        if ws === nothing
            vs = f(xs...)
            if !primalcalled
                value = primal_value(lowest(ab), vs, f, xs)
                primalcalled = true
            end
            return value, nothing
        end
        pb = pullback_function(lowest(ab), (_xs...,) -> begin
            vs = f(_xs...)
            if !primalcalled
                value = primal_value(lowest(ab), vs, f, xs)
                primalcalled = true
            end
            return vs
        end, xs...)(ws)
        return value, pb
    end
end

struct LazyDerivative{B, F, X}
    backend::B
    f::F
    xs::X
end

function Base.:*(d::LazyDerivative, y)
    return derivative(d.backend, d.f, d.xs...) * y
end

function Base.:*(y, d::LazyDerivative)
    return y * derivative(d.backend, d.f, d.xs...)
end

function Base.:*(d::LazyDerivative, y::Union{Number,Tuple})
    if y isa Tuple && d.xs isa Tuple
        @assert length(y) == length(d.xs)
    end
    return derivative(d.backend, d.f, d.xs...) .* y
end

function Base.:*(y::Union{Number,Tuple}, d::LazyDerivative)
    if y isa Tuple && d.xs isa Tuple
        @assert length(y) == length(d.xs)
    end
    return y .* derivative(d.backend, d.f, d.xs...)
end

function Base.:*(d::LazyDerivative, y::AbstractArray)
    return map((d)-> d*y, derivative(d.backend, d.f, d.xs...))
end

function Base.:*(y::AbstractArray, d::LazyDerivative)
    return map((d)-> y*d, derivative(d.backend, d.f, d.xs...))
end

struct LazyGradient{B, F, X}
    backend::B
    f::F
    xs::X
end
Base.:*(d::LazyGradient, y) = gradient(d.backend, d.f, d.xs...) * y
Base.:*(y, d::LazyGradient) = y * gradient(d.backend, d.f, d.xs...)

function Base.:*(d::LazyGradient, y::Union{Number,Tuple})
    if y isa Tuple && d.xs isa Tuple
        @assert length(y) == length(d.xs)
    end
    if d.xs isa Tuple
        return gradient(d.backend, d.f, d.xs...) .* y
    else
        return gradient(d.backend, d.f, d.xs) .* y
    end
end

function Base.:*(y::Union{Number,Tuple}, d::LazyGradient)
    if y isa Tuple && d.xs isa Tuple
        @assert length(y) == length(d.xs)
    end
    if d.xs isa Tuple
        return y .* gradient(d.backend, d.f, d.xs...)
    else
        return y .* gradient(d.backend, d.f, d.xs)
    end
end

struct LazyJacobian{B, F, X}
    backend::B
    f::F
    xs::X
end

function Base.:*(d::LazyJacobian, ys)
    if !(ys isa Tuple)
        ys = (ys, )
    end
    if d.xs isa Tuple
        vjp = pushforward_function(d.backend, d.f, d.xs...)(ys)
    else
        vjp = pushforward_function(d.backend, d.f, d.xs)(ys)
    end
    if vjp isa Tuple
        return vjp
    else
        return (vjp,)
    end
end

function Base.:*(ys, d::LazyJacobian)
    if ys isa Tuple
        ya = adjoint.(ys)
    else
        ya = adjoint(ys)
    end
    if d.xs isa Tuple
        return pullback_function(d.backend, d.f, d.xs...)(ya)
    else
        return pullback_function(d.backend, d.f, d.xs)(ya)
    end
end

function Base.:*(d::LazyJacobian, ys::Number)
    if d.xs isa Tuple
        return jacobian(d.backend, d.f, d.xs...) .* ys
    else
        return jacobian(d.backend, d.f, d.xs) .* ys
    end
end

function Base.:*(ys::Number, d::LazyJacobian)
    if d.xs isa Tuple
        return jacobian(d.backend, d.f, d.xs...) .* ys
    else
        return jacobian(d.backend, d.f, d.xs) .* ys
    end
end

struct LazyHessian{B, F, X}
    backend::B
    f::F
    xs::X
end

function Base.:*(d::LazyHessian, ys)
    if !(ys isa Tuple)
        ys = (ys, )
    end

    if d.xs isa Tuple
        res =  pushforward_function(
            second_lowest(d.backend),
            (xs...,) -> gradient(lowest(d.backend), d.f, xs...)[1], d.xs...,)(ys)  # [1] because gradient returns a tuple
    else
        res =  pushforward_function(
            second_lowest(d.backend),
            (xs,) -> gradient(lowest(d.backend), d.f, xs)[1],d.xs,)(ys)  # gradient returns a tuple
    end
    if res isa Tuple
        return res
    else
        return (res,)
    end
end

function Base.:*(ys, d::LazyHessian)
    if ys isa Tuple
        ya = adjoint.(ys)
    else
        ya = adjoint(ys)
    end
    if d.xs isa Tuple
        return pullback_function(
            second_lowest(d.backend),
            (xs...,) -> gradient(lowest(d.backend), d.f, xs...),
            d.xs...,
            )(ya)
    else
        return pullback_function(
            second_lowest(d.backend),
            (xs,) -> gradient(lowest(d.backend), d.f, xs)[1],
            d.xs,
            )(ya)
    end
end

function Base.:*(d::LazyHessian, ys::Number)
    if d.xs isa Tuple
        return hessian(d.backend, d.f, d.xs...).*ys
    else
        return hessian(d.backend, d.f, d.xs).*ys
    end
end

function Base.:*(ys::Number, d::LazyHessian)
    if d.xs isa Tuple
        return ys.*hessian(d.backend, d.f, d.xs...)
    else
        return ys.*hessian(d.backend, d.f, d.xs)
    end
end

function lazy_derivative(ab::AbstractBackend, f, xs::Number...)
    return LazyDerivative(ab, f, xs)
end
function lazy_gradient(ab::AbstractBackend, f, xs...)
    return LazyGradient(ab, f, xs)
end
function lazy_hessian(ab::AbstractBackend, f, xs...)
    return LazyHessian(ab, f, xs)
end
function lazy_jacobian(ab::AbstractBackend, f, xs...)
    return LazyJacobian(ab, f, xs)
end

struct D{B, F}
    backend::B
    f::F
end
D(b::AbstractBackend, d::D) = H(HigherOrderBackend((b, d.b)), d.f)
D(d::D) = H(HigherOrderBackend((d.backend, d.backend)), d.f)
function (d::D)(xs...; lazy = true)
    if lazy
        return lazy_jacobian(d.ab, d.f, xs...)
    else
        return jacobian(d.ab, d.f, xs...)
    end
end

struct H{B, F}
    backend::B
    f::F
end
function (h::H)(xs...; lazy = true)
    if lazy
        return lazy_hessian(h.ab, h.f, xs...)
    else
        return hessian(h.ab, h.f, xs...)
    end
end

macro primitive(expr)
    fdef = ExprTools.splitdef(expr)
    name = fdef[:name]
    if name == :pushforward_function
        return define_pushforward_function_and_friends(fdef) |> esc
    elseif name == :pullback_function
        return define_pullback_function_and_friends(fdef) |> esc
    elseif name == :jacobian
        return define_jacobian_and_friends(fdef) |> esc
    elseif name == :primal_value
        return define_primal_value(fdef) |> esc
    else
        throw("Unsupported AD primitive.")
    end
end

function define_pushforward_function_and_friends(fdef)
    fdef[:name] = :(AbstractDifferentiation.pushforward_function)
    args = fdef[:args]
    funcs = quote
        $(ExprTools.combinedef(fdef))
        function AbstractDifferentiation.jacobian($(args...),)
            identity_like = AbstractDifferentiation.identity_matrix_like($(args[3:end]...),)
            pff = AbstractDifferentiation.pushforward_function($(args...),)
            if eltype(identity_like) <: Tuple{Vararg{Union{AbstractMatrix, Number}}}
                return map(identity_like) do identity_like_i
                    if VERSION < v"1.3"
                        return reduce(hcat, map(AbstractDifferentiation._eachcol.(identity_like_i)...) do (cols...)
                            pff(cols)
                        end)
                    else
                    return mapreduce(hcat, AbstractDifferentiation._eachcol.(identity_like_i)...) do (cols...)
                        pff(cols)
                    end
                end
                end
            elseif eltype(identity_like) <: AbstractMatrix
                # needed for the computation of the Hessian and Jacobian
                ret = hcat.(mapslices(identity_like[1], dims=1) do cols
                    # cols loop over basis states
                    pf = pff((cols,))
                    if typeof(pf) <: AbstractVector
                        # to make the hcat. work / get correct matrix-like, non-flat output dimension
                        return (pf, )
                    else
                        return pf
                    end
                end ...)
                return ret isa Tuple ? ret : (ret,)

            else
                return pff(identity_like)
            end
        end
    end
    return funcs
end

function define_pullback_function_and_friends(fdef)
    fdef[:name] = :(AbstractDifferentiation.pullback_function)
    args = fdef[:args]
    funcs = quote
        $(ExprTools.combinedef(fdef))
        function AbstractDifferentiation.jacobian($(args...),)
            value_and_pbf = AbstractDifferentiation.value_and_pullback_function($(args...),)
            value, _ = value_and_pbf(nothing)
            identity_like = AbstractDifferentiation.identity_matrix_like(value)
            if eltype(identity_like) <: Tuple{Vararg{AbstractMatrix}}
                return map(identity_like) do identity_like_i
                    if VERSION < v"1.3"
                        return reduce(vcat, map(AbstractDifferentiation._eachcol.(identity_like_i)...) do (cols...)
                            value_and_pbf(cols)[2]'
                        end)
                    else
                    return mapreduce(vcat, AbstractDifferentiation._eachcol.(identity_like_i)...) do (cols...)
                        value_and_pbf(cols)[2]'
                        end
                    end
                end
            elseif eltype(identity_like) <: AbstractMatrix
                # needed for Hessian computation:
                # value is a (grad,). Then, identity_like is a (matrix,).
                # cols loops over columns of the matrix
                return vcat.(mapslices(identity_like[1], dims=1) do cols
                    adjoint.(value_and_pbf((cols,))[2])
                end ...)
            else
                return adjoint.(value_and_pbf(identity_like)[2])
            end
        end
    end
    return funcs
end

_eachcol(a::Number) = (a,)
_eachcol(a) = eachcol(a)

function define_jacobian_and_friends(fdef)
    fdef[:name] = :(AbstractDifferentiation.jacobian)
    return ExprTools.combinedef(fdef)
end

function define_primal_value(fdef)
    fdef[:name] = :(AbstractDifferentiation.primal_value)
    return ExprTools.combinedef(fdef)
end

function identity_matrix_like(x)
    throw("The function `identity_matrix_like` is not defined for the type $(typeof(x)).")
end
function identity_matrix_like(x::AbstractVector)
    return (Matrix{eltype(x)}(I, length(x), length(x)),)
end
function identity_matrix_like(x::Number)
    return (one(x),)
end
identity_matrix_like(x::Tuple) = identity_matrix_like(x...)
@generated function identity_matrix_like(x...)
    expr = :(())
    for i in 1:length(x)
        push!(expr.args, :(()))
        for j in 1:i-1
            push!(expr.args[i].args, :((zero_matrix_like(x[$j])[1])))
        end
        push!(expr.args[i].args, :((identity_matrix_like(x[$i]))[1]))
        for j in i+1:length(x)
            push!(expr.args[i].args, :(zero_matrix_like(x[$j])[1]))
        end
    end
    return expr
end

zero_matrix_like(x::Tuple) = zero_matrix_like(x...)
zero_matrix_like(x...) = map(zero_matrix_like, x)
zero_matrix_like(x::AbstractVector) = (zero(similar(x, length(x), length(x))),)
zero_matrix_like(x::Number) = (zero(x),)
function zero_matrix_like(x)
    throw("The function `zero_matrix_like` is not defined for the type $(typeof(x)).")
end

function __init__()
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")
end

end


text:
"""
DatingEstimate is a simple datatype to store the results of the coaltime function.
coaltime estimates date range that may be considered a 95% confidence range in
which the true coalsescence time between two sequences lies.
The type stores the 5%, 50%, and 95% values for the range.
"""
struct DatingEstimate
    lower::Float64
    middle::Float64
    upper::Float64
end

DatingEstimate() = DatingEstimate(.0, .0, .0)

lower(x::DatingEstimate) = x.lower
middle(x::DatingEstimate) = x.middle
upper(x::DatingEstimate) = x.upper

lower(x::Vector{DatingEstimate}) = lower.(x)
middle(x::Vector{DatingEstimate}) = middle.(x)
upper(x::Vector{DatingEstimate}) = upper.(x)

function Base.show(io::IO, de::DatingEstimate)
    println(io, "Coalescence time estimate:\n5%: $(de.lower), 95%: $(de.upper)")
end

@inline function binomzero(p0::Float64, N::Int, B::Int)
    f(p::Float64) = cdf(Binomial(N, p), B) - p0
    return fzero(f, 0, 1)
end

"""
    estimate_time(N::Int, K::Int, µ::Float64)

Compute the coalescence time between two sequences by modelling the process of
mutation accumulation between two sequences as a bernoulli process.

`N` is the length of the two aligned sequences.
`K` is the estimated number of mutations.
`µ` is the assumed mutation rate.
"""
@inline function estimate_time(N::Int, K::Int, µ::Float64)
    if !(N ≥ K ≥ 0)
        throw(ArgumentError("Condition must be satisfied: N ≥ K ≥ 0"))
    end
    div = 2µ
    ninetyfive = ceil(binomzero(0.05, N, K) / div)
    fifty = ceil(binomzero(0.5, N, K) / div)
    five = ceil(binomzero(0.95, N, K) / div)
    return DatingEstimate(five, fifty, ninetyfive)
end

"""
    estimate_time(N::Int, K::Int, µ::Float64)

Compute the coalescence time between two sequences by modelling the process of
mutation accumulation between two sequences as a bernoulli process.

`N` is the length of the two aligned sequences.
`d` is the a value of estimated genetic distance.
`µ` is the assumed mutation rate.
"""
@inline function estimate_time(N::Int, d::Float64, µ::Float64)
    if !(1.0 ≥ d ≥ 0.0)
        throw(ArgumentError("Genetic distance `d` must be a value between 1 and 0"))
    end
    K = convert(Int, ceil(d * N))
    return estimate_time(N, K, µ)
end

@inline function estimate_time(m::MutationCount, μ::Float64)
    return estimate_time(m.n_sites, m.n_mutations, μ)
end

@inline estimate_time(V::Vector{MutationCount}, μ::Float64) = estimate_time.(V, μ)
@inline estimate_time(V::Vector{Vector{MutationCount}}, μ::Float64) = estimate_time.(V, μ)

@inline function estimate_time(M::PairwiseListMatrix{T,false}, μ::Float64) where {T}
    V = estimate_time(getlist(M), μ)
    return PairwiseListMatrix(V, false, eltype(V)())
end

@inline function estimate_time(M::N, μ::Float64) where {N<:NamedArray}
    V = estimate_time(M.array, μ)
    return setlabels(V, names(M, 1))
end


text:
module simulator

export matrixcore, fresnelltools, analyticalmaterials, cmtgratings, utilities, materials

include("utilities.jl")
include("analyticalmaterials.jl")
include("ema.jl")
include("materials.jl")
include("matrixcore.jl")
include("fresnelltools.jl")
include("targetfigures.jl")

end # simulator


text:
import LinearAlgebra.BLAS: gemv, gemv!, gemm!, trsm!, axpy!, ger!
import LinearAlgebra: cholesky, Cholesky
import Base: getproperty

Base.@deprecate chol(X) cholesky(X).U

const AM = AbstractMatrix
const UT = UpperTriangular

@explicit_intercepts(
    Cholesky,
    Tuple{AbstractMatrix{<:∇Scalar}, Union{Char, Symbol}, Integer},
    [true, false, false],
)
function ∇(
    ::Type{Cholesky},
    ::Type{Arg{1}},
    p,
    C::Cholesky,
    X̄::Union{UpperTriangular, LowerTriangular},
    X::Union{UpperTriangular, LowerTriangular},
    uplo::Union{Char, Symbol},
    info::Integer,
)
    # We aren't doing any actual computation if we've constructed a Cholesky object
    # directly, so just pass through this call and return the sensitivies
    return X̄
end
function ∇(
    ::Type{Cholesky},
    ::Type{Arg{1}},
    p,
    C::Cholesky,
    X̄::Tangent{<:Cholesky},
    X::Union{UpperTriangular, LowerTriangular},
    uplo::Union{Char, Symbol},
    info::Integer,
)
    return getproperty(X̄, Symbol(uplo))
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   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.

"""
subtract.jl defines methods for the -(x, y) function
"""

# --- Exports

import Base.:(-)

# --- Method definitions

-(M::AbstractMultivector, N::AbstractMultivector) = M + -N

# ------ Specializations involving a Pseudoscalar instance

# B::Pseudoscalar, C::Pseudoscalar
-(B::Pseudoscalar, C::Pseudoscalar) = Pseudoscalar(B, value=value(B) - value(C))

# ------ Specializations involving an AbstractScalar instance

# B::AbstractScalar, C::One
# B::One, C::AbstractScalar
-(B::AbstractScalar, C::One) = Scalar{typeof(value(B))}(value(B) - 1)
-(B::One, C::AbstractScalar) = Scalar{typeof(value(C))}(1 - value(C))

# B::AbstractScalar, C::Zero
# B::Zero, C::AbstractScalar
-(B::AbstractScalar, C::Zero) = B
-(B::Zero, C::AbstractScalar) = -C

# B::AbstractScalar, x::Real
# x::Real, B::AbstractScalar
-(B::AbstractScalar, x::Real) = Scalar{typeof(value(B))}(value(B) - x)
-(x::Real, B::AbstractScalar) = Scalar{typeof(value(B))}(x - value(B))

# ------ Specializations involving a Scalar instance

# B::Scalar, C::Scalar
-(B::Scalar, C::Scalar) = Scalar{typeof(value(B))}(value(B) - value(C))

# ------ Specializations involving a One instance

# B::One, C::One
-(B::One, C::One) = Zero{typeof(value(B))}()

# B::One, C::Zero
# B::Zero, C::One
-(B::One, C::Zero) = B
-(B::Zero, C::One) = -C

# ------ Specializations involving a Zero instance

# B::Zero, C::Zero
-(B::Zero, C::Zero) = B

# B::Zero, v::Vector
# v::Vector, B::Zero
-(B::Zero, v::Vector{<:Real}) = -v
-(v::Vector{<:Real}, B::Zero) = v


text:
##############
# InnerLabel #
##############
immutable InnerProduct{N} <: Number
    b::StateLabel{N}
    k::StateLabel{N}
end

blabel(i::InnerProduct) = i.b
klabel(i::InnerProduct) = i.k

Base.repr(i::InnerProduct) = brstr(blabel(i))*ktstr(klabel(i))[2:end]
Base.show(io::IO, i::InnerProduct) = print(io, repr(i))

Base.(:(==))(a::InnerProduct, b::InnerProduct) = blabel(a) == blabel(b) && klabel(a) == klabel(b)
Base.(:(==))(::InnerProduct, ::Number) = false
Base.(:(==))(::Number, ::InnerProduct) = false

Base.hash(i::InnerProduct) = hash(blabel(i), hash(klabel(i)))
Base.hash(i::InnerProduct, h::Uint64) = hash(hash(i), h)

Base.conj(i::InnerProduct) = InnerProduct(klabel(i), blabel(i))

##############
# inner_rule #
##############

# we can cheat here to avoid tempory StateLabel construction
# to evaluate KroneckerDelta inner products
eval_inner_rule(p::KroneckerDelta, b, k) = inner_rule(p, b, k)
eval_inner_rule(p::KroneckerDelta, b::StateLabel, k::StateLabel) = inner_rule(p, b, k)

eval_inner_rule(p::AbstractInner, b::StateLabel, k::StateLabel) = inner_rule(p, b, k)
eval_inner_rule(p::AbstractInner, b, k) = inner_rule(p, StateLabel(b), StateLabel(k))

inner_rule{P<:AbstractInner}(::P, b, k) = error("inner_rule(::$P, ::StateLabel, ::StateLabel) must be defined to evaluate inner products with type $P")
inner_rule(::UndefinedInner, b, k) = InnerExpr(InnerProduct(b, k))
inner_rule(::KroneckerDelta, b, k) = b == k ? 1 : 0

inner_rettype{P<:AbstractInner}(::P) = first(Base.return_types(inner_rule, (P, StateLabel, StateLabel)))
inner_rettype(::UndefinedInner) = InnerExpr
inner_rettype(::KroneckerDelta) = Int64

# very common operation for state/operator inner products
inner_mul(v,c,prodtype,b,k) = v * c * eval_inner_rule(prodtype, b, k)

##############
# InnerExpr #
##############
# A InnerExpr is a type that wraps arthimetic expressions
# performed with InnerExprs. The allows storage and
# delayed evaluation of expressions. For example, this
# expression:
#
#   (< a | b >^2 + < c | d >^2 - 3.13+im) / 2
#
# is representable as a InnerExpr.

immutable InnerExpr <: Number
    ex::Expr
end

InnerExpr(iex::InnerExpr) = InnerExpr(iex.ex)
InnerExpr{N<:Number}(n::N) = convert(InnerExpr, n)

Base.convert(::Type{InnerExpr}, iex::InnerExpr) = iex
Base.convert{N<:Number}(::Type{InnerExpr}, n::N) = InnerExpr(Expr(:call, +, n))

Base.(:(==))(a::InnerExpr, b::InnerExpr) = a.ex == b.ex
same_num(a::Number, b::Number) = a == b
same_num(a::InnerExpr, b::InnerExpr) = a == b
same_num(iex::InnerExpr, i::InnerProduct) = iex == InnerExpr(i)
same_num(i::InnerProduct, iex::InnerExpr) = iex == i
same_num(iex::InnerExpr, n::Number) = iex == InnerExpr(n)
same_num(n::Number, iex::InnerExpr) = iex == n

Base.hash(iex::InnerExpr) = hash(iex.ex)
Base.hash(iex::InnerExpr, h::Uint64) = hash(hash(iex), h)

Base.one(::InnerExpr) = InnerExpr(1)
Base.zero(::InnerExpr) = InnerExpr(0)

Base.promote_rule{N<:Number}(::Type{InnerExpr}, ::Type{N}) = InnerExpr

Base.length(iex::InnerExpr) = length(iex.ex.args)
Base.getindex(iex::InnerExpr, i) = iex.ex.args[i]

##############
# inner_eval #
##############
inner_eval(f, iex::InnerExpr) = eval(inner_reduce!(f, copy(iex.ex)))
inner_eval(f, c) = c

inner_reduce!(f, iex::InnerExpr) = inner_reduce!(f, copy(iex.ex))
inner_reduce!(f::Function, i::InnerProduct) = f(blabel(i), klabel(i))
inner_reduce!(p::AbstractInner, i::InnerProduct) = inner_rule(p, blabel(i), klabel(i))
inner_reduce!(f, n) = n

function inner_reduce!(f, ex::Expr)
    for i=1:length(ex.args)
        ex.args[i] = inner_reduce!(f, ex.args[i])
    end
    return ex
end

######################
# Printing Functions #
######################
Base.show(io::IO, iex::InnerExpr) = print(io, repr(iex.ex)[2:end])

##################
# Exponentiation #
##################
function iexpr_exp(a, b)
    if same_num(b, 1)
        return InnerExpr(a)
    elseif same_num(a, 0)
        return InnerExpr(0)
    elseif same_num(b, 0)
        return InnerExpr(1)
    else
        return InnerExpr(:($(s)^$(n)))
    end
end

Base.(:^)(a::InnerExpr, b::Integer) = iexpr_exp(a, b)
Base.(:^)(a::InnerExpr, b::Rational) = iexpr_exp(a, b)
Base.(:^)(a::InnerExpr, b::InnerExpr) = iexpr_exp(a, b)
Base.(:^)(a::InnerExpr, b::Number) = iexpr_exp(a, b)
Base.(:^)(a::MathConst{:e}, b::InnerExpr) = iexpr_exp(a, b)
Base.(:^)(a::Number, b::InnerExpr) = iexpr_exp(a, b)

Base.exp(iex::InnerExpr) = InnerExpr(:(exp($(iex))))
Base.exp2(iex::InnerExpr) = InnerExpr(:(exp2($(iex))))

Base.sqrt(iex::InnerExpr) = InnerExpr(:(sqrt($(iex))))

Base.log(iex::InnerExpr) = length(iex)==2 && iex[1]==:exp ? iex[2] : InnerExpr(:(log($(iex))))
Base.log(a::MathConst{:e}, b::InnerExpr) = InnerExpr(:(log($(a),$(b))))
Base.log(a::InnerExpr, b::InnerExpr) = InnerExpr(:(log($(a),$(b))))
Base.log(a::InnerExpr, b::Number) = InnerExpr(:(log($(a),$(b))))
Base.log(a::Number, b::InnerExpr) = InnerExpr(:(log($(a),$(b))))
Base.log2(iex::InnerExpr) = length(iex)==2 && iex[1]==:exp2 ? iex[2] : InnerExpr(:(log2($(iex))))

##################
# Multiplication #
##################
function iexpr_mul(a,b)
    if same_num(b, 1)
        return InnerExpr(a)
    elseif same_num(a, 1)
        return InnerExpr(b)
    elseif same_num(b, 0) || same_num(a, 0)
        return InnerExpr(0)
    else
        return InnerExpr(:($(a)*$(b)))
    end
end

Base.(:*)(a::InnerExpr, b::InnerExpr) =iexpr_mul(a,b)
Base.(:*)(a::Bool, b::InnerExpr) = iexpr_mul(a,b)
Base.(:*)(a::InnerExpr, b::Bool) = iexpr_mul(a,b)
Base.(:*)(a::InnerExpr, b::Number) = iexpr_mul(a,b)
Base.(:*)(a::Number, b::InnerExpr) = iexpr_mul(a,b)

############
# Division #
############
function iexpr_div(a, b)
    if same_num(a, b)
        return InnerExpr(1)
    elseif same_num(a, 0)
        return InnerExpr(0)
    elseif same_num(b, 0)
        return InnerExpr(Inf)
    elseif same_num(b, 1)
        return InnerExpr(a)
    else
        return InnerExpr(:($(a)/$(b)))
    end
end

Base.(:/)(a::InnerExpr, b::InnerExpr) = iexpr_div(a, b)
Base.(:/)(a::InnerExpr, b::Complex) = iexpr_div(a, b)
Base.(:/)(a::InnerExpr, b::Number) = iexpr_div(a, b)
Base.(:/)(a::Number, b::InnerExpr) = iexpr_div(a, b)

############
# Addition #
############
function iexpr_add(a, b)
    if same_num(a, 0)
        return InnerExpr(b)
    elseif same_num(b, 0)
        return InnerExpr(a)
    else
        return InnerExpr(:($(a)+$(b)))
    end
end

Base.(:+)(a::InnerExpr, b::InnerExpr) = iexpr_add(a, b)
Base.(:+)(a::InnerExpr, b::Number) = iexpr_add(a, b)
Base.(:+)(a::Number, b::InnerExpr) = iexpr_add(a, b)

###############
# Subtraction #
###############
Base.(:-)(iex::InnerExpr) = length(iex)==2 && iex[1]==:- ? InnerExpr(iex[2]) :  InnerExpr(:(-$(iex)))

function iexpr_subtract(a, b)
    if same_num(a, b)
        return InnerExpr(0)
    elseif same_num(a, 0)
        return InnerExpr(-b)
    elseif same_num(b, 0)
        return InnerExpr(a)
    else
        return InnerExpr(:($(a)-$(b)))
    end
end

Base.(:-)(a::InnerExpr, b::InnerExpr) = iexpr_subtract(a, b)
Base.(:-)(a::InnerExpr, b::Number) = iexpr_subtract(a, b)
Base.(:-)(a::Number, b::InnerExpr) = iexpr_subtract(a, b)

##################
# Absolute Value #
##################
Base.abs(iex::InnerExpr) = length(iex)==2 && iex[1]==:abs ? iex :  InnerExpr(:(abs($(iex))))
Base.abs2(iex::InnerExpr) = InnerExpr(:(abs2($iex)))

#####################
# Complex Conjugate #
#####################
Base.conj(iex::InnerExpr) = length(iex)==2 && iex[1]==:conj ? InnerExpr(iex[2]) :  InnerExpr(:(conj($(iex))))
Base.ctranspose(iex::InnerExpr) = conj(iex)

##########################
# Elementwise Operations #
##########################
for op=(:*,:-,:+,:/,:^)
    elop = symbol(string(:.) * string(op))
    @eval begin
        ($elop)(a::InnerExpr, b::InnerExpr) = ($op)(a,b)
        ($elop)(a::InnerExpr, b::Number) = ($op)(a,b)
        ($elop)(a::Number, b::InnerExpr) = ($op)(a,b)
    end
end

######################
# Printing Functions #
######################
function Base.repr(iex::InnerExpr)
    if iex[1]==+
        if length(iex)==2
            return repr(iex[2])
        else
            return repr(Expr(:call, iex.ex[2:end]...))[3:end-1]
        end
    else
        return repr(iex.ex)[2:end]
    end
end

Base.show(io::IO, iex::InnerExpr) = print(io, repr(iex))

export InnerExpr,
    inner_eval,
    inner_rule,
    inner_rettype


text:
@testset "misc" begin
    @testset "bypass_graph" begin
        N = 4
        E = 5
        adj = [0 1 1 1;
            1 0 1 0;
            1 1 0 1;
            1 0 1 0]

        nf = rand(3, N)
        ef = rand(5, E)
        gf = rand(7)

        fg = FeaturedGraph(adj, nf=nf, ef=ef, gf=gf)
        layer = bypass_graph(x -> x .+ 1.,
                                x -> x .+ 2.,
                                x -> x .+ 3.)
        fg_ = layer(fg)
        @test GraphSignals.adjacency_matrix(fg_) == adj
        @test node_feature(fg_) == nf .+ 1.
        @test edge_feature(fg_) == ef .+ 2.
        @test global_feature(fg_) == gf .+ 3.
    end
end


text:
using RecipesBase
using ColorTypes, ColorSchemes
using TiledIteration

@recipe function atmimshow(cc::Cnsolt{T,2}; cscheme=ColorSchemes.gray, rangescale=(-0.5,0.5), atomscale=10, coordinate=:cartesian) where {T}
    nch = nchannels(cc)
    ord = orders(cc)
    df =  decimations(cc)

    layout :=  (2,nch)
    size -->  100 .* df .* (ord .+ 1) .* (nch, 2)

    atmsreim = map(analysiskernels(cc)) do fp
        apply_colorscheme(Val{coordinate}, fp, cscheme, rangescale, atomscale)
    end

    atmsup = [ [ f[1] for f in atmsreim ]..., [ f[2] for f in atmsreim ]... ]

    for idx = 1:2*nch
        @series begin
            subplot := idx
            axis    := false
            grid    := false

            aspect_ratio := :equal
            atmsup[idx]
        end
    end
end

@recipe function atmimshow(cc::Cnsolt{T,2}, p::Integer; cscheme=ColorSchemes.gray, rangescale=(-0.5,0.5), atomscale=10, coordinate=:cartesian) where {T}
    nch = nchannels(cc)
    ord = orders(cc)
    df =  decimations(cc)

    layout :=  (2,1)
    size -->  100 .* df .* (ord .+ 1) .* (1,2)

    afs = analysiskernels(cc)

    atmsup = apply_colorscheme(Val{coordinate}, afs[p], cscheme, rangescale, atomscale)

    for idx = 1:2
        @series begin
            subplot := idx
            axis    := false
            grid    := false

            aspect_ratio := :equal
            atmsup[idx]
        end
    end
end

# @recipe function atmimshow(cc::Cnsolt{T,2}, ::Val{:onecolumn}; kwargs...)
#     atmimshow(cc; kwargs...)
# end

@recipe function atmimshow(cc::Rnsolt{T,2}; cscheme=ColorSchemes.gray, rangescale=(-0.5,0.5), atomscale=10) where {T}
    mxP = maximum(cc.nChannels)
    nch = cc.nChannels
    ord = orders(cc)
    df =  decimations(cc)
    difch = nch[2]-nch[1]

    layout     :=  (2,mxP)
    size -->  100 .* df .* (ord .+ 1) .* (mxP, 2)

    afs = analysiskernels(cc)
    dummyimg = fill(-Inf, size(afs[1]))
    afssym = [ afs[1:nch[1]]; fill(dummyimg, max(difch, 0)) ]
    afsasym = [ afs[(nch[1]+1):end]; fill(dummyimg, max(-difch, 0)) ]

    afsup, afslw = map((afssym, afsasym,)) do fs
        map(f->apply_colorscheme(f, cscheme, rangescale, atomscale), fs)
    end

    for idx = 1:mxP
        @series begin
            subplot := idx
            axis    := false
            grid    := false
            aspect_ratio := :equal

            afsup[idx]
        end

        @series begin
            subplot := idx+mxP
            axis    := false
            grid    := false
            aspect_ratio := :equal

            afslw[idx]
        end
    end
end

@recipe function atmimshow(cc::Rnsolt{T,2}, p::Integer; cscheme=ColorSchemes.gray, rangescale=(-0.5,0.5), atomscale=10) where {T}
    ord = orders(cc)
    df =  decimations(cc)

    axis    := false
    grid    := false

    aspect_ratio := :equal
    size -->  100 .* df .* (ord .+ 1)

    afs = analysiskernels(cc)

    apply_colorscheme(afs[p], cscheme, rangescale, atomscale)
end

@recipe function atmimshow(cc::Rnsolt{T,2}, ::Val{:onecolumn}; cscheme=ColorSchemes.gray, rangescale=(-0.5,0.5), atomscale=10) where {T}
    mxP = maximum(cc.nChannels)
    ord = orders(cc)
    df =  decimations(cc)
    nch = nchannels(cc)

    layout     :=  (1, nch)
    size -->  200 .* df .* (ord .+ 1) .* (nch, 1)

    afsout = map(f->apply_colorscheme(f, cscheme, rangescale, atomscale), analysiskernels(cc))

    for idx = 1:nch
        @series begin
            subplot := idx
            axis    := false
            grid    := false
            aspect_ratio := :equal

            afsout[idx]
        end
    end
end

apply_colorscheme(x::AbstractArray, args...; kwargs...) = apply_colorscheme(Val{:cartesian}, x, args...; kwargs...)
apply_colorscheme(::Type{Val{:cartesian}}, x::AbstractArray{T}, cscheme::AbstractVector{C}, rangescale::Tuple{R,R}, atomscale::Integer) where {T<:Complex,C<:Colorant,R<:Real} = map(t->apply_colorscheme(t, cscheme, rangescale, atomscale), reim(x))
function apply_colorscheme(::Type{Val{:cartesian}}, x::AbstractArray{T}, cscheme::AbstractVector{C}, rangescale::Tuple{R,R}, atomscale::Integer) where {T<:Real,C<:Colorant,R<:Real}
    resize_by_nn(get(cscheme, x, rangescale), atomscale)
end

function apply_colorscheme(::Type{Val{:polar}}, x::AbstractArray{T}, cscheme::AbstractVector{C}, rangescale::Tuple{R,R}, atomscale::Integer) where {T,C<:Colorant,R<:Real}
    mxv = norm(rangescale)
    rx = get(cscheme, abs.(x), (0.0, mxv))
    ax = RGB.(HSV.(x .|> angle .|> rad2deg .|> wrapdeg, 1.0, abs.(x) / mxv))
    resize_by_nn.((rx, ax), atomscale)
end

function resize_by_nn(x::AbstractArray{T,D}, scale::Integer) where {T,D}
    output = similar(x, T, size(x) .* scale)

    for (idx, tile) in enumerate(TileIterator(axes(output), (fill(scale,D)...,)))
        output[tile...] .= x[idx]
    end
    output
end

# wrapdeg(x::Real) = ifelse(x < 0, x+360, x)
wrapdeg(x::Real) = mod(x, 360)

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
using NoiseContrastiveEstimation
using Documenter
ENV["GKSwstype"] = "100"
@show joinpath(@__DIR__, "src")

makedocs(;
    modules = [NoiseContrastiveEstimation],
    authors = "Francesco Alemanno <francescoalemanno710@gmail.com> and contributors",
    repo = "https://github.com/francescoalemanno/NoiseContrastiveEstimation.jl/blob/{commit}{path}#L{line}",
    sitename = "NoiseContrastiveEstimation.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://francescoalemanno.github.io/NoiseContrastiveEstimation.jl",
        assets = String[],
    ),
    pages = ["Basic Usage" => "index.md", "Reference" => "reference.md"],
)

deploydocs(;
    repo = "github.com/francescoalemanno/NoiseContrastiveEstimation.jl",
    push_preview = true,
)