text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@everywhere module TestDecisionFunctions

using StochasticPrograms
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MA = MOI.MutableArithmetics

function test_vectorization(x, fx, y, fy, z, fz, w, fw)
    g = VectorAffineDecisionFunction(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], [z, w])),
            [3, 1, 4],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([3, 6], [x, y])),
            zeros(Int, 3),
        ),
    )
    @testset "vectorize" begin
        g1 = AffineDecisionFunction(
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2, w)], 3),
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(6, y)], 0),
        )
        g2 = AffineDecisionFunction(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Int}[], 1),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Int}[], 0),
        )
        g3 = AffineDecisionFunction(
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(5, z)], 4),
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 0),
        )
        @test g ≈ MOIU.vectorize([g1, g2, g3])
        vov = MOIU.vectorize(SingleDecision[])
        @test MOI.output_dimension(vov) == 0
        @test vov isa VectorOfDecisions
        aff = MOIU.vectorize(AffineDecisionFunction{Int}[])
        @test MOI.output_dimension(aff) == 0
        @test aff isa VectorAffineDecisionFunction{Int}
        # TODO
        # quad = MOIU.vectorize(QuadraticDecisionFunction{Int}[])
        # @test MOI.output_dimension(quad) == 0
        # @test quad isa VectorQuadraticDecisionFunction{Int}
    end
    @testset "operate vcat" begin
        d = VectorOfDecisions([x, y])
        @testset "Decision with $T" for T in [Int, Float64]
            @test VectorOfDecisions == MOIU.promote_operation(
                vcat,
                T,
                typeof(fx),
                typeof(d),
                typeof(fy),
            )
            vod = MOIU.operate(vcat, T, fx, d, fy)
            @test vod.decisions == [x, x, y, y]
            @test VectorOfDecisions == MOIU.promote_operation(
                vcat,
                T,
                typeof(d),
                typeof(fx),
                typeof(fy),
            )
            vod = MOIU.operate(vcat, T, d, fx, fy)
            @test vod.decisions == [x, y, x, y]
            @test VectorOfDecisions == MOIU.promote_operation(
                vcat,
                T,
                typeof(fx),
                typeof(fy),
                typeof(d),
            )
            vod = MOIU.operate(vcat, T, fx, fy, d)
            @test vod.decisions == [x, y, x, y]
        end
        v = MOI.VectorOfVariables([z, w])
        f = AffineDecisionFunction(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2, 4], [z, w]), 5),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5, 3], [x, y]), 0),
        )
        g = VectorAffineDecisionFunction(
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.(
                    [3, 1],
                    MOI.ScalarAffineTerm.([5, 2], [z, w]),
                ),
                [3, 1, 4],
            ),
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.(
                    [3, 1],
                    MOI.ScalarAffineTerm.([3, 6], [x, y]),
                ),
                zeros(Int, 3),
            ),
        )
        @testset "Affine" begin
            @test MOIU.promote_operation(
                vcat,
                Int,
                typeof(fx),
                typeof(f),
                typeof(d),
                typeof(fz),
                typeof(v),
                Int,
                typeof(g),
                typeof(fy),
                Int,
                typeof(fw),
            ) == VectorAffineDecisionFunction{Int}
            F = MOIU.operate(vcat, Int, fx, f, d, fz, v, 3, g, fy, -4, fw)
            expected_variable_terms =
                MOI.VectorAffineTerm.(
                    [2, 2, 5, 6, 7, 11, 9, 14],
                    MOI.ScalarAffineTerm.(
                        [2, 4, 1, 1, 1, 5, 2, 1],
                        [z, w, z, z, w, z, w, w],
                    ),
                )
            expected_decision_terms =
                MOI.VectorAffineTerm.(
                    [1, 2, 2, 3, 4, 11, 9, 12],
                    MOI.ScalarAffineTerm.(
                        [1, 5, 3, 1, 1, 3, 6, 1],
                        [x, x, y, x, y, x, y, y],
                    ),
                )
            expected_constants = [0, 5, 0, 0, 0, 0, 0, 3, 3, 1, 4, 0, -4, 0]
            @test F.variable_part.terms == expected_variable_terms
            @test F.decision_part.terms == expected_decision_terms
            @test F.variable_part.constants == expected_constants
        end
    end
end

function test_eval_variables(x, fx, y, fy, z, fz, w, fw)
    vals = Dict(x => 3, y => 1, z => 5, w => 1)
    @test MOI.output_dimension(fx) == 1
    @test MOIU.eval_variables(vi -> vals[vi], fx) ≈ 3
    @test MOIU.eval_variables(vi -> vals[vi], fx) ≈ 3
    fxy = VectorOfDecisions([x, y])
    @test MOI.output_dimension(fxy) == 2
    @test MOIU.eval_variables(vi -> vals[vi], fxy) ≈ [3, 1]
    @test MOIU.eval_variables(vi -> vals[vi], fxy) ≈ [3, 1]
    f = AffineDecisionFunction(
        MOI.ScalarAffineFunction(
            [
                MOI.ScalarAffineTerm(3.0, z),
                MOI.ScalarAffineTerm(2.0, w),
            ],
            2.0,
        ),
        MOI.ScalarAffineFunction(
            [
                MOI.ScalarAffineTerm(1.0, x),
                MOI.ScalarAffineTerm(2.0, y),
            ],
            0.0,
        ),
    )
    @test MOI.output_dimension(f) == 1
    @test MOIU.eval_variables(vi -> vals[vi], f) ≈ 24
    @test MOIU.eval_variables(vi -> vals[vi], f) ≈ 24
    f = VectorAffineDecisionFunction(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 2],
                MOI.ScalarAffineTerm.([1.0, 2.0], [z, w]),
            ),
            [-3.0, 2.0],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 2],
                MOI.ScalarAffineTerm.([3.0, 2.0], [x, y]),
            ),
            zeros(Float64, 2),
        ),
    )
    @test MOI.output_dimension(f) == 2
    @test MOIU.eval_variables(vi -> vals[vi], f) ≈ [6, 11]
    @test MOIU.eval_variables(vi -> vals[vi], f) ≈ [6, 11]
    f = QuadraticDecisionFunction(
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm.(1.0, [z, w]),
            MOI.ScalarQuadraticTerm.(1.0, [z, w], [z, w]),
            -12.0,
        ),
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm.(1.0, [x, y]),
            MOI.ScalarQuadraticTerm.(1.0, [x, x, y], [x, y, y]),
            0.0,
        ),
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm{Float64}[],
            MOI.ScalarQuadraticTerm.(1.0, [z, w], [x, y]),
            0.0,
        ),
    )
    @test MOI.output_dimension(f) == 1
    @test MOIU.eval_variables(vi -> vals[vi], f) ≈ 35
    @test MOIU.eval_variables(vi -> vals[vi], f) ≈ 35
    # TODO
    # f = VectorQuadraticDecisionFunction(
    #     MOI.VectorQuadraticFunction(
    #         MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.(1.0, [x, y])),
    #         MOI.VectorQuadraticTerm.(
    #             [1, 2, 2],
    #             MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]),
    #         ),
    #         [-3.0, -2.0],
    #     ),
    #     MOI.VectorQuadraticFunction(
    #         MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.(1.0, [x, y])),
    #         MOI.VectorQuadraticTerm.(
    #             [1, 2, 2],
    #             MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]),
    #         ),
    #         [-3.0, -2.0],
    #     ),
    #     MOI.VectorQuadraticFunction(
    #         MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.(1.0, [x, y])),
    #         MOI.VectorQuadraticTerm.(
    #             [1, 2, 2],
    #             MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]),
    #         ),
    #         [-3.0, -2.0],
    #     ),
    # )
    # @test MOI.output_dimension(fvq) == 2
    # @test MOIU.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1]
    # @test MOIU.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1]
end

function test_substitute_variables(x, fx, y, fy, z, fz, w, fw)
    # We do tests twice to make sure the function is not modified
    subs = Dict(x => 2.0fy + 1.0, y => 1.0fy, z => -1.0fw, w => 1.0fy + 1.0fz)
    vals = Dict(x => 3.0, y => 1.0, z => 5.0, w => 0.0)
    subs_vals = Dict(x => 3.0, y => 1.0, z => 0.0, w => 6.0)
    # Affine
    f = fx + 2.0fy + 3.0fz + 2.0
    f_subbed = -3.0fw + 4.0fy + 3.0
    @test MOIU.eval_variables(vi -> subs_vals[vi], f) ==
        MOIU.eval_variables(vi -> vals[vi], f_subbed)
    @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
    @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
    # VectorAffine
    f = MOIU.operate(vcat, Float64, 3.0fz - 3.0, fx + 2.0fy + 2.0)
    f_subbed = MOIU.operate(vcat, Float64, -3.0fw - 3.0, 4.0fy + 3.0)
    @test MOIU.eval_variables(vi -> subs_vals[vi], f) ==
        MOIU.eval_variables(vi -> vals[vi], f_subbed)
    @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
    @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
    # Quadratic
    f = 1.0fx + 1.0fy + 1.0fx * fz + 1.0fw * fz + 1.0fw * fy + 2.0fw * fw - 3.0
    f_subbed =
        3.0fy - 1.0fw + 3.0fy * fy + 2.0fz * fz - 3.0fy * fw - 1.0fz * fw +
        5.0fy * fz - 2.0
    @test MOIU.eval_variables(vi -> subs_vals[vi], f) ==
        MOIU.eval_variables(vi -> vals[vi], f_subbed)
    @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
    @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
    # TODO
    # f = MOIU.operate(
    #     vcat,
    #     Float64,
    #     1.0fy + 1.0fx * fz - 3.0,
    #     1.0fx + 1.0fw * fz + 1.0fw * fy - 2.0,
    # )
    # f_subbed = MOIU.operate(
    #     vcat,
    #     Float64,
    #     1.0fy - fw - 2.0fy * fw - 3.0,
    #     2.0fy + 1.0fy * fy - 1.0fw * fy - 1.0fw * fz + 1.0fy * fz - 1.0,
    # )
    # @test MOIU.eval_variables(vi -> subs_vals[vi], f) ==
    #       MOIU.eval_variables(vi -> vals[vi], f_subbed)
    # @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
    # @test MOIU.substitute_variables(vi -> subs[vi], f) ≈ f_subbed
end

function test_map_indices(x, fx, y, fy, z, fz, w, fw)
    f = QuadraticDecisionFunction(
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm.(1.0, [z, w]),
            MOI.ScalarQuadraticTerm.(1.0, [z, w], [z, w]),
            -12.0,
        ),
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm.(1.0, [x, y]),
            MOI.ScalarQuadraticTerm.(1.0, [x, x, y], [x, y, y]),
            0.0,
        ),
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm{Float64}[],
            MOI.ScalarQuadraticTerm.(1.0, [z, w], [x, y]),
            0.0,
        ),
    )
    index_map = Dict(x => y, y => x, w => w, z => w)
    g = MOIU.map_indices(index_map, f)
    @test g.variable_part.affine_terms ==
        MOI.ScalarAffineTerm.(1.0, [w, w])
    @test g.variable_part.quadratic_terms ==
        MOI.ScalarQuadraticTerm.(1.0, [w, w], [w, w])
    @test g.decision_part.affine_terms ==
        MOI.ScalarAffineTerm.(1.0, [y, x])
    @test g.decision_part.quadratic_terms ==
        MOI.ScalarQuadraticTerm.(1.0, [y, y, x], [y, x, x])
    @test isempty(g.cross_terms.affine_terms)
    @test g.cross_terms.quadratic_terms ==
        MOI.ScalarQuadraticTerm.(1.0, [w, w], [y, x])
    @test MOI.constant(g) == -12.0
end

function test_VectorOfDecisions_iteration(x, fx, y, fy, z, fz, w, fw)
    f = VectorOfDecisions([x, y, x])
    it = MOIU.eachscalar(f)
    @test length(it) == 3
    @test eltype(it) == SingleDecision
    @test collect(it) == [
        SingleDecision(x),
        SingleDecision(y),
        SingleDecision(x),
    ]
    @test it[2] == SingleDecision(y)
    @test it[end] == SingleDecision(x)
end

function test_VectorAffineDecisionFunction_iteration(x, fx, y, fy, z, fz, w, fw)
    f = VectorAffineDecisionFunction(
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [3, 2, 1],
                MOI.ScalarAffineTerm.(
                    [2, 9, 1],
                    [z, z, w],
                ),
            ),
            [2, 7, 5],
        ),
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [2, 2, 2, 3, 1, 2],
                MOI.ScalarAffineTerm.(
                    [1, 7, 3, 6, 4, 1],
                    [x, y, y, x, x, y],
                ),
            ),
            zeros(Int, 3),
        ),
    )
    it = MOIU.eachscalar(f)
    @test length(it) == 3
    @test eltype(it) == AffineDecisionFunction{Int}
    g = it[2]
    @test g isa AffineDecisionFunction{Int}
    @test g.variable_part.terms ==
        MOI.ScalarAffineTerm.([9], [z])
    @test g.decision_part.terms ==
        MOI.ScalarAffineTerm.([1, 7, 3, 1], [x, y, y, y])
    @test MOI.constant(g) == 7
    g = it[1]
    @test g isa AffineDecisionFunction{Int}
    @test g.variable_part.terms ==
        MOI.ScalarAffineTerm.([1], [w])
    @test g.decision_part.terms ==
        MOI.ScalarAffineTerm.([4], [x])
    @test MOI.constant(g) == 2
    g = it[end]
    @test g isa AffineDecisionFunction{Int}
    @test g.variable_part.terms ==
        MOI.ScalarAffineTerm.([2], [z])
    @test g.decision_part.terms ==
        MOI.ScalarAffineTerm.([6], [x])
    @test MOI.constant(g) == 5
    h = it[[2, 1]]
    @test h isa VectorAffineDecisionFunction{Int}
    @test sort(h.variable_part.terms, by = t -> t.output_index) ==
        MOI.VectorAffineTerm.(
            [1, 2],
            MOI.ScalarAffineTerm.([9, 1], [z, w]),
        )
    @test sort(h.decision_part.terms, by = t -> t.output_index) ==
        MOI.VectorAffineTerm.(
            [1, 1, 1, 1, 2],
            MOI.ScalarAffineTerm.([1, 7, 3, 1, 4], [x, y, y, y, x]),
        )
    @test MOIU.constant_vector(h) == [7, 2]
    F = MOIU.operate(vcat, Int, it[[1, 2]], it[3])
    @test F isa VectorAffineDecisionFunction{Int}
    @test sort(F.variable_part.terms, by = t -> t.output_index) ==
        MOI.VectorAffineTerm.(
            [1, 2, 3],
            MOI.ScalarAffineTerm.(
                [1, 9, 2],
                [w, z, z],
            ),
        )
    @test sort(F.decision_part.terms, by = t -> t.output_index) ==
        MOI.VectorAffineTerm.(
            [1, 2, 2, 2, 2, 3],
            MOI.ScalarAffineTerm.(
                [4, 1, 7, 3, 1, 6],
                [x, x, y, y, y, x],
            ),
        )
    @test MOIU.constant_vector(F) == MOIU.constant_vector(f)
end

function test_scalar_operations(x, fx, y, fy, z, fz, w, fw)
    @testset "Decision" begin
        f = SingleDecision(MOI.VariableIndex(0))
        g = SingleDecision(MOI.VariableIndex(1))
        @testset "one" begin
            @test !isone(f)
            @test !isone(g)
        end
        @testset "zero" begin
            @test !iszero(f)
            @test !iszero(g)
            @test f + 1 ≈ 1 + f
            @test (f + 1.0) - 1.0 ≈ (2.0f) / 2.0
            @test (f - 1.0) + 1.0 ≈ (2.0f) / 2.0
            @test (1.0 + f) - 1.0 ≈ (f * 2.0) / 2.0
            @test 1.0 - (1.0 - f) ≈ (f / 2.0) * 2.0
        end
    end
    @testset "Affine" begin
        @testset "zero" begin
            f = @inferred MOIU.zero(AffineDecisionFunction{Float64})
            @test iszero(f)
            @test MOIU.isapprox_zero(f, 1e-16)
        end
        @testset "promote_operation" begin
            @test MOIU.promote_operation(-, Int, SingleDecision) ==
                AffineDecisionFunction{Int}
            @test MOIU.promote_operation(
                -,
                Int,
                AffineDecisionFunction{Int},
            ) == AffineDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Float64,
                SingleDecision,
                SingleDecision,
            ) == AffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Float64,
                MOI.SingleVariable,
                SingleDecision,
            ) == AffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Float64,
                SingleDecision,
                MOI.SingleVariable,
            ) == AffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Float64,
                AffineDecisionFunction{Float64},
                Float64,
            ) == AffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Int,
                AffineDecisionFunction{Int},
                AffineDecisionFunction{Int},
            ) == AffineDecisionFunction{Int}
        end
        @testset "Comparison" begin
            @test MOIU.operate(
                +,
                Float64,
                SingleDecision(x),
                SingleDecision(y),
            ) + 1.0 ≈ AffineDecisionFunction(
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm{Float64}[],
                    1.0,
                ),
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm.([1, 1e-7, 1], [x, y, y]),
                    0.0,
                ),
            ) atol = 1e-6
            f1 = AffineDecisionFunction(
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm(1.0, z), MOI.ScalarAffineTerm(1e-7, w)],
                    1.0,
                ),
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(1e-7, y)],
                    0.0,
                ),
            )
            f2 = AffineDecisionFunction(
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, z)], 1.0),
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            )
            @test f1 ≈ f2 atol = 1e-6
            fdiff = f1 - f2
            @testset "With iszero" begin
                @test !iszero(fdiff)
                @test iszero(f1 - f1)
                @test iszero(f2 - f2)
            end
            @testset "With tolerance" begin
                MOIU.canonicalize!(fdiff)
                @test !MOIU.isapprox_zero(fdiff, 1e-8)
                @test MOIU.isapprox_zero(fdiff, 1e-6)
            end
        end
        @testset "canonical" begin
            f = MOIU.canonical(
                AffineDecisionFunction(
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([1, 1, 3, -2, -3], [w, z, z, w, z]),
                        5,
                    ),
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([2, 1, 3, -2, -3], [y, x, y, x, y]),
                        0,
                    ),
                )
            )
            @test MOI.output_dimension(f) == 1
            @test f.variable_part.terms == MOI.ScalarAffineTerm.([1, -1], [z, w])
            @test f.decision_part.terms == MOI.ScalarAffineTerm.([-1, 2], [x, y])
            @test MOI.constant(f) == 5
        end
        f = AffineDecisionFunction(
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([0.5, 0.5], [z, w]),
                0.5,
            ),
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 0.5], [x, y]),
                0.0,
            ),
        )
        @testset "convert" begin
            @test_throws InexactError convert(SingleDecision, f)
            @test_throws InexactError MOIU.convert_approx(SingleDecision, f)
            @test MOIU.convert_approx(SingleDecision, f, tol = 0.5) ==
                SingleDecision(x)
            @test convert(typeof(f), f) === f
            quad_f = QuadraticDecisionFunction(
                MOI.ScalarQuadraticFunction(
                    f.variable_part.terms,
                    MOI.ScalarQuadraticTerm{Float64}[],
                    f.variable_part.constant,
                ),
                MOI.ScalarQuadraticFunction(
                    f.decision_part.terms,
                    MOI.ScalarQuadraticTerm{Float64}[],
                    0.0,
                ),
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarAffineTerm{Float64}[],
                    MOI.ScalarQuadraticTerm{Float64}[],
                    0.0,
                ),
            )
            @test convert(QuadraticDecisionFunction{Float64}, f) ≈ quad_f
            for g in [
                convert(
                    AffineDecisionFunction{Float64},
                    SingleDecision(x),
                ),
                convert(
                    AffineDecisionFunction{Float64},
                    1SingleDecision(x),
                ),
            ]
                @test g isa AffineDecisionFunction{Float64}
                @test convert(SingleDecision, g) == SingleDecision(x)
                @test MOIU.convert_approx(SingleDecision, g) ==
                    SingleDecision(x)
            end
        end
        @testset "operate with Float64 coefficient type" begin
            f = AffineDecisionFunction(
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm.([1.0, 4.0], [z, w]),
                    5.0,
                ),
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm.([1.0, 4.0], [x, y]),
                    0.0,
                ),
            )
            @test f ≈ 2.0f / 2.0
        end
        @testset "operate with Int coefficient type" begin
            f = MOIU.canonical(
                AffineDecisionFunction(
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.(
                            [1, 1, 2],
                            [w, w, z],
                        ),
                        2,
                    ),
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.(
                            [3, 2, -3],
                            [y, x, x],
                        ),
                        0,
                    ),
                ) + AffineDecisionFunction(
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([-2, -2], [z, w]),
                        3,
                    ),
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([-1, 3, 2], [y, x, y]),
                        0,
                    ),
                )
            )
            @test f === +f
            @test f ≈
                SingleDecision(x) + AffineDecisionFunction(
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        5,
                    ),
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([1, 4], [x, y]),
                        0,
                    ),
                )
            @test f ≈ f * 1
            @test f ≈
                AffineDecisionFunction(
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        2,
                    ),
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([1, 2], [x, y]),
                        0,
                    ),
                ) * 2 + 1
            @test f ≈
                SingleDecision(x) - AffineDecisionFunction(
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        -5,
                    ),
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([-1, -4], [x, y]),
                        -5,
                    ),
                )
            @test f ≈
                AffineDecisionFunction(
                    MOI.ScalarAffineFunction(
                        [MOI.ScalarAffineTerm(1, w)],
                        5,
                    ),
                    MOI.ScalarAffineFunction(
                        MOI.ScalarAffineTerm.([3, 4], [x, y]),
                        5,
                    ),
                ) - SingleDecision(x) - MOI.SingleVariable(w)
        end
        @testset "modification" begin
            f = MOIU.modify_function(f, MOI.ScalarConstantChange(6))
            @test MOI.constant(f) == 6
            g = deepcopy(f)
            @test g ≈ f
            f = MOIU.modify_function(f, MOI.ScalarCoefficientChange(z, 3))
            f = MOIU.modify_function(f, StochasticPrograms.DecisionCoefficientChange(y, 3))
            @test !(g ≈ f)
            @test g.variable_part.terms == MOI.ScalarAffineTerm[]
            @test g.decision_part.terms == MOI.ScalarAffineTerm.([2, 4], [x, y])
            @test f.variable_part.terms == [MOI.ScalarAffineTerm(3, z)]
            @test f.decision_part.terms == MOI.ScalarAffineTerm.([2, 3], [x, y])
            f = MOIU.modify_function(f, StochasticPrograms.DecisionCoefficientChange(x, 0))
            @test f.decision_part.terms == MOI.ScalarAffineTerm.([3], [y])
        end
    end
    @testset "Quadratic" begin
        @testset "zero" begin
            f = @inferred MOIU.zero(QuadraticDecisionFunction{Float64})
            @test MOIU.isapprox_zero(f, 1e-16)
        end
        @testset "promote_operation" begin
            @test MOIU.promote_operation(
                -,
                Int,
                QuadraticDecisionFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Int,
                QuadraticDecisionFunction{Int},
                QuadraticDecisionFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Int,
                QuadraticDecisionFunction{Int},
                MOI.ScalarAffineFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Int,
                MOI.ScalarAffineFunction{Int},
                QuadraticDecisionFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Int,
                QuadraticDecisionFunction{Int},
                MOI.ScalarQuadraticFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Int,
                MOI.ScalarQuadraticFunction{Int},
                QuadraticDecisionFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Int,
                QuadraticDecisionFunction{Int},
                AffineDecisionFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Int,
                AffineDecisionFunction{Int},
                QuadraticDecisionFunction{Int},
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                *,
                Int,
                SingleDecision,
                SingleDecision,
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                *,
                Int,
                MOI.SingleVariable,
                SingleDecision,
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                *,
                Int,
                SingleDecision,
                MOI.SingleVariable,
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                *,
                Float64,
                SingleDecision,
                MOI.ScalarAffineFunction{Float64},
            ) == QuadraticDecisionFunction{Float64}
            @test MOIU.promote_operation(
                *,
                Int,
                MOI.ScalarAffineFunction{Int},
                SingleDecision,
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                *,
                Float64,
                MOI.SingleVariable,
                AffineDecisionFunction{Float64},
            ) == QuadraticDecisionFunction{Float64}
            @test MOIU.promote_operation(
                *,
                Int,
                AffineDecisionFunction{Int},
                MOI.SingleVariable,
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                *,
                Float64,
                SingleDecision,
                AffineDecisionFunction{Float64},
            ) == QuadraticDecisionFunction{Float64}
            @test MOIU.promote_operation(
                *,
                Int,
                AffineDecisionFunction{Int},
                SingleDecision,
            ) == QuadraticDecisionFunction{Int}
            @test MOIU.promote_operation(
                *,
                Float64,
                MOI.ScalarAffineFunction{Float64},
                AffineDecisionFunction{Float64},
            ) == QuadraticDecisionFunction{Float64}
            @test MOIU.promote_operation(
                *,
                Float64,
                AffineDecisionFunction{Float64},
                MOI.ScalarAffineFunction{Float64},
            ) == QuadraticDecisionFunction{Float64}
            @test MOIU.promote_operation(
                *,
                Float64,
                AffineDecisionFunction{Float64},
                AffineDecisionFunction{Float64},
            ) == QuadraticDecisionFunction{Float64}
            @test MOIU.promote_operation(
                /,
                Float64,
                QuadraticDecisionFunction{Float64},
                Float64,
            ) == QuadraticDecisionFunction{Float64}
        end
        f = 7 + 3fx + 1fx * fx + 2fy * fy + 3fx * fy + 1fz * fz + 1fx * fz
        MOIU.canonicalize!(f)
        @test MOI.output_dimension(f) == 1
        @testset "Comparison" begin
            @testset "With iszero" begin
                @test !iszero(f)
                @test iszero(0 * f)
                @test iszero(f - f)
            end
            @testset "With tolerance" begin
                @test !MOIU.isapprox_zero(f, 1e-8)
                @test MOIU.isapprox_zero(0 * f, 1e-8)
                g = 1.0fx * fy - (1 + 1e-6) * fy * fx + 1.0fx * fz - (1 + 1e-6) * fz * fx
                MOIU.canonicalize!(g)
                @test MOIU.isapprox_zero(g, 1e-5)
                @test !MOIU.isapprox_zero(g, 1e-7)
            end
        end
        @testset "convert" begin
            @test_throws InexactError convert(SingleDecision, f)
            @test_throws InexactError convert(AffineDecisionFunction{Int}, f)
            g = convert(QuadraticDecisionFunction{Float64}, fx)
            @test convert(SingleDecision, g) == fx
        end
        @testset "Power" begin
            @testset "Affine" begin
                aff = 1fx + 2 + fy + 1fz
                @test isone(@inferred aff^0)
                @test convert(typeof(f), aff) ≈ @inferred aff^1
                @test aff * aff ≈ @inferred aff^2
                err =
                    ArgumentError("Cannot take $(typeof(aff)) to the power 3.")
                @test_throws err aff^3
            end
            @testset "Quadratic" begin
                @test isone(@inferred f^0)
                @test f ≈ @inferred f^1
                err = ArgumentError("Cannot take $(typeof(f)) to the power 2.")
                @test_throws err f^2
            end
        end
        @testset "operate" begin
            @testset "No zero affine term" begin
                qd(f) = convert(QuadraticDecisionFunction{Int}, f)
                for fzfw in [qd(1fz * fw), qd(fw * 1fz)]
                    @test isempty(fzfw.variable_part.affine_terms)
                    @test length(fzfw.variable_part.quadratic_terms) == 1
                    @test isempty(fzfw.decision_part.affine_terms)
                    @test isempty(fzfw.decision_part.quadratic_terms)
                    @test fzfw.variable_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, z, w) ||
                        fzfw.variable_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, w, z)
                end
                for fzfz in [qd(1fz * fz), qd(fz * 1fz)]
                    @test isempty(fzfz.variable_part.affine_terms)
                    @test length(fzfz.variable_part.quadratic_terms) == 1
                    @test isempty(fzfz.decision_part.affine_terms)
                    @test isempty(fzfz.decision_part.quadratic_terms)
                    @test fzfz.variable_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(2, z, z)
                end
                for fxfy in [1fx * fy, fx * 1fy]
                    @test isempty(fxfy.variable_part.affine_terms)
                    @test isempty(fxfy.variable_part.quadratic_terms)
                    @test isempty(fxfy.decision_part.affine_terms)
                    @test length(fxfy.decision_part.quadratic_terms) == 1
                    @test fxfy.decision_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, x, y) ||
                        fxfy.decision_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, y, x)
                end
                for fxfx in [1fx * fx, fx * 1fx]
                    @test isempty(fxfx.variable_part.affine_terms)
                    @test isempty(fxfx.variable_part.quadratic_terms)
                    @test isempty(fxfx.decision_part.affine_terms)
                    @test length(fxfx.decision_part.quadratic_terms) == 1
                    @test fxfx.decision_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(2, x, x)
                end
                for fxfy in [1fx * fy, fx * 1fy]
                    @test isempty(fxfy.decision_part.affine_terms)
                    @test length(fxfy.decision_part.quadratic_terms) == 1
                    @test fxfy.decision_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, x, y) ||
                        fxfy.decision_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, y, x)
                end
                for fxfx in [1fx * fx, fx * 1fx]
                    @test isempty(fxfx.decision_part.affine_terms)
                    @test length(fxfx.decision_part.quadratic_terms) == 1
                    @test fxfx.decision_part.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(2, x, x)
                end
                for fxfz in [1fx * fz, fz * 1fx]
                    @test isempty(fxfz.variable_part.affine_terms)
                    @test isempty(fxfz.variable_part.quadratic_terms)
                    @test isempty(fxfz.decision_part.affine_terms)
                    @test isempty(fxfz.decision_part.quadratic_terms)
                    @test isempty(fxfz.cross_terms.affine_terms)
                    @test length(fxfz.cross_terms.quadratic_terms) == 1
                    @test fxfz.cross_terms.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, x, z) ||
                        fxfz.cross_terms.quadratic_terms[1] ==
                        MOI.ScalarQuadraticTerm(1, z, x)
                end
            end
            @testset "operate!" begin
                q = 1.0fx + 1.0fy + (1.0fx) * fy + (1.0fx) * fz + (1.0fw) * fz
                @test q ≈ 1.0fx + 1.0fy + (1.0fw) * fz + (1.0fx) * fz + (1.0fx) * fy
                # This calls
                aff = 1.0fx + 1.0fy
                # which tries to mutate `aff`, gets a quadratic expression
                # and mutate it with the remaining term
                @test MOIU.operate!(
                    +,
                    Float64,
                    aff,
                    (1.0fx) * fz,
                    (1.0fw) * fz,
                ) ≈ q + (1.0fx) * fy
            end
            f = 7 + 3fx + 2fx * fx + 2fy * fy + 3fx * fy + 3fx * fz + 2fz * fz
            @test f ≈ 7 + (fx + 2fy) * (1fx + fy) + 3fx + (fx + 2fz) * (1fx + fz)
            @test f ≈ -(-7 - 3fx) + (fx + 2fy) * (1fx + fy) + (fx + 2fz) * (1fx + fz)
            @test f ≈ -((fx + 2fy) * (MOIU.operate(-, Int, fx) - fy)) + 3fx + 7 - ((fx + 2fz) * (MOIU.operate(-, Int, fx) - fz))
            @test f ≈
                7 + MOIU.operate(*, Int, fx, fx) + 3fx * (fy + 1) + 2fy * fy +
                MOIU.operate(*, Int, fx, fx) + 3fx * fz + 2fz * fz
            @test f ≈
                (fx + 2) * (fx + 1) + (fy + 1) * (2fy + 3fx) + (fx + 2) * (fx + 1) + (fz + 1) * (2fz + 3fx) + (3 - 9fx - 2fy - 2fz)
            @test f ≈ begin
                QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([4], [z], [z]),
                        4,
                    ),
                    MOI.ScalarQuadraticFunction(
                        [MOI.ScalarAffineTerm(3, x)],
                        MOI.ScalarQuadraticTerm.([4], [x], [x]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm{Int}[],
                        0,
                    ),
                ) + QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm{Int}[],
                        3,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([4, 3], [y, x], [y, y]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([3], [x], [z]),
                        0,
                    ),
                )
            end
            @test f ≈ begin
                QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([4], [z], [z]),
                        10,
                    ),
                    MOI.ScalarQuadraticFunction(
                        [MOI.ScalarAffineTerm(3, x)],
                        MOI.ScalarQuadraticTerm.([2], [x], [x]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm{Int}[],
                        0,
                    ),
                ) - QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm{Int}[],
                        3,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([-2, -4, -3], [x, y, x], [x, y, y]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([-3], [x], [z]),
                        0,
                    ),
                )
            end
            @test f ≈ begin
                AffineDecisionFunction(
                    MOI.ScalarAffineFunction{Int}(
                        MOI.ScalarAffineTerm{Int}[],
                        5,
                    ),
                    MOI.ScalarAffineFunction{Int}(
                        [MOI.ScalarAffineTerm(3, x)],
                        0,
                    ),
                ) + QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        [MOI.ScalarAffineTerm(2, z)],
                        MOI.ScalarQuadraticTerm.([4], [z], [z]),
                        2,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([4, 4, 3], [x, y, x], [x, y, y]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([3], [x], [z]),
                        0,
                    ),
                ) - MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm(2, z)],
                    0,
                )
            end
            @test f ≈ begin
                AffineDecisionFunction(
                    MOI.ScalarAffineFunction{Int}(
                        MOI.ScalarAffineTerm{Int}[],
                        5,
                    ),
                    MOI.ScalarAffineFunction{Int}(
                        [MOI.ScalarAffineTerm(3, x)],
                        0,
                    ),
                ) - QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        [MOI.ScalarAffineTerm(2, z)],
                        MOI.ScalarQuadraticTerm.([-4], [z], [z]),
                        -2,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([-4, -4, -3], [x, y, x], [x, y, y]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([3], [x], [z]),
                        0,
                    ),
                ) + MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm(2, z)],
                    0,
                )
            end
            @test f ≈ begin
                QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([4], [z], [z]),
                        2,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([4, 4, 3], [x, y, x], [x, y, y]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([3], [x], [z]),
                        0,
                    ),
                ) + AffineDecisionFunction(
                    MOI.ScalarAffineFunction{Int}(
                        MOI.ScalarAffineTerm{Int}[],
                        5,
                    ),
                    MOI.ScalarAffineFunction{Int}(
                        [MOI.ScalarAffineTerm(3, x)],
                        0,
                    ),
                )
            end
            @test f ≈ begin
                QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        [MOI.ScalarAffineTerm(2, z)],
                        MOI.ScalarQuadraticTerm.([4], [z], [z]),
                        9,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([4, 4, 3], [x, y, x], [x, y, y]),
                        0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Int}[],
                        MOI.ScalarQuadraticTerm.([3], [x], [z]),
                        0,
                    ),
                ) - AffineDecisionFunction(
                    MOI.ScalarAffineFunction{Int}(
                        MOI.ScalarAffineTerm{Int}[],
                        5,
                    ),
                    MOI.ScalarAffineFunction{Int}(
                        [MOI.ScalarAffineTerm(-2, x)],
                        0,
                    ),
                ) - 2 * fz + 1 * fx + 3
            end
            @test f ≈ begin
                2.0 * QuadraticDecisionFunction(
                    MOI.ScalarQuadraticFunction(
                        [MOI.ScalarAffineTerm(2., z)],
                        MOI.ScalarQuadraticTerm.([4.], [z], [z]),
                        7.0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        [MOI.ScalarAffineTerm(3., x)],
                        MOI.ScalarQuadraticTerm.([4., 4., 3.], [x, y, x], [x, y, y]),
                        0.0,
                    ),
                    MOI.ScalarQuadraticFunction(
                        MOI.ScalarAffineTerm{Float64}[],
                        MOI.ScalarQuadraticTerm.([3.], [x], [z]),
                        0.0,
                    ),
                ) / 2.0
            end
            @testset "modification" begin
                f = MOIU.modify_function(f, MOI.ScalarConstantChange(9))
                @test MOI.constant(f) == 9
                f = MOIU.modify_function(f, MOI.ScalarCoefficientChange(z, 0))
                f = MOIU.modify_function(f, StochasticPrograms.DecisionCoefficientChange(y, 0))
                @test isempty(f.variable_part.affine_terms)
                @test f.decision_part.affine_terms == MOI.ScalarAffineTerm.([3], [x])
                g = deepcopy(f)
                @test f ≈ g
                f = MOIU.modify_function(f, MOI.ScalarCoefficientChange(z, 2))
                f = MOIU.modify_function(f, StochasticPrograms.DecisionCoefficientChange(y, 2))
                @test !(f ≈ g)
                @test isempty(g.variable_part.affine_terms)
                @test f.variable_part.affine_terms == MOI.ScalarAffineTerm.([2], [z])
                @test g.decision_part.affine_terms == MOI.ScalarAffineTerm.([3], [x])
                @test f.decision_part.affine_terms == MOI.ScalarAffineTerm.([3, 2], [x, y])
            end
        end
    end
end

function test_vector_operations(x, fx, y, fy, z, fz, w, fw)
    @testset "Affine" begin
        @testset "promote_operation" begin
            @test MOIU.promote_operation(-, Int, VectorOfDecisions) ==
                VectorAffineDecisionFunction{Int}
            @test MOIU.promote_operation(
                -,
                Int,
                VectorAffineDecisionFunction{Int},
            ) == VectorAffineDecisionFunction{Int}
            @test MOIU.promote_operation(
                +,
                Float64,
                VectorOfDecisions,
                VectorOfDecisions,
            ) == VectorAffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Float64,
                MOI.VectorOfVariables,
                VectorOfDecisions,
            ) == VectorAffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Float64,
                VectorOfDecisions,
                MOI.VectorOfVariables,
            ) == VectorAffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Float64,
                VectorAffineDecisionFunction{Float64},
                Vector{Float64},
            ) == VectorAffineDecisionFunction{Float64}
            @test MOIU.promote_operation(
                +,
                Int,
                VectorAffineDecisionFunction{Int},
                VectorAffineDecisionFunction{Int},
            ) == VectorAffineDecisionFunction{Int}
        end
        @testset "Comparison" begin
            @test MOIU.operate(
                +,
                Float64,
                MOI.VectorOfVariables([z, w]),
                VectorOfDecisions([x, y]),
            ) + ones(2) ≈ VectorAffineDecisionFunction(
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1., 1.], [z, w])),
                    [1., 1.],
                ),
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1, 2, 2], MOI.ScalarAffineTerm.([1., 1e-7, 1.], [x, y, y])),
                    zeros(Float64, 2),
                ),
            ) atol = 1e-6
            f1 = VectorAffineDecisionFunction(
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1.0, 1e-7], [z, w])),
                    [1., 1.],
                ),
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1., 1e-7], [x, y])),
                    zeros(Float64, 2),
                ),
            )
            f2 = VectorAffineDecisionFunction(
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1], MOI.ScalarAffineTerm.([1.], [z])),
                    [1., 1.],
                ),
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1], MOI.ScalarAffineTerm.([1.], [x])),
                    zeros(Float64, 2),
                ),
            )
            @test f1 ≈ f2 atol = 1e-6
        end
        @testset "canonical" begin
            f = MOIU.canonical(
                VectorAffineDecisionFunction(
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([2, 1, 1, 2, 1], MOI.ScalarAffineTerm.([1, 1, 3, -2, -3], [w, z, z, w, z])),
                        [5, 2],
                    ),
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([2, 1, 2, 1, 2], MOI.ScalarAffineTerm.([2, 1, 3, -2, -3], [y, x, y, x, y])),
                        zeros(Int, 2),
                    ),
                )
            )
            @test MOI.output_dimension(f) == 2
            @test f.variable_part.terms == MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1, -1], [z, w]))
            @test f.decision_part.terms == MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([-1, 2], [x, y]))
            @test MOI.constant(f) == [5, 2]
        end
        f = VectorAffineDecisionFunction(
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([0.5, 0.5], [z, w])),
                [0.5, 0.5],
            ),
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1., 1.], [x, y])),
                zeros(Float64, 2),
            ),
        )
        @testset "convert" begin
            @test_throws InexactError MOIU.convert_approx(VectorOfDecisions, f)
            @test MOIU.convert_approx(VectorOfDecisions, f, tol = 0.5) ==
                VectorOfDecisions([x, y])
        end
        @testset "operate with Float64 coefficient type" begin
            f = VectorAffineDecisionFunction(
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1.0, 4.0], [z, w])),
                    [5.0, 2.0],
                ),
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1., 4.], [x, y])),
                    zeros(Float64, 2),
                ),
            )
            @test f ≈ 2.0f / 2.0
        end
        @testset "operate with Int coefficient type" begin
            f = MOIU.canonical(
                VectorAffineDecisionFunction(
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 1, 2], MOI.ScalarAffineTerm.([1, 1, 2], [w, w, z])),
                        [5, 0],
                    ),
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 1, 2], MOI.ScalarAffineTerm.([2, -3, 3], [x, x, y])),
                        zeros(Int, 2),
                    ),
                ) + VectorAffineDecisionFunction(
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([-2, -2], [w, z])),
                        [0, 2],
                    ),
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([2, 1, 2], MOI.ScalarAffineTerm.([-1, 3, 2], [y, x, y])),
                        zeros(Int, 2),
                    ),
                )
            )
            @test f === +f
            @test f ≈
                VectorOfDecisions([x, x]) + VectorAffineDecisionFunction(
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm{Int}[],
                        [5, 2],
                    ),
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 2, 2], MOI.ScalarAffineTerm.([1, 4, -1], [x, y, x])),
                        zeros(Int, 2),
                    ),
                )
            @test f ≈ f * 1
            @test f ≈
                VectorAffineDecisionFunction(
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm{Int}[],
                        [0, 1],
                    ),
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1, 2], [x, y])),
                        zeros(Int, 2),
                    ),
                ) * 2 + [5, 0]
            @test f ≈
                VectorOfDecisions([x, y]) - VectorAffineDecisionFunction(
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm{Int}[],
                        [-5, -2],
                    ),
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([-1, -3], [x, y])),
                        zeros(Int, 2),
                    ),
                )
            @test f ≈ MOIU.canonical(
                VectorAffineDecisionFunction(
                    MOI.VectorAffineFunction(
                        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1, w))],
                        [5, 2],
                    ),
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([3, 4], [x, y])),
                        zeros(Int, 2),
                    ),
                ) - MOIU.operate(vcat, Int, 1SingleDecision(x) + MOI.SingleVariable(w), 0)
            )
        end
        @testset "modification" begin
            f = MOIU.modify_function(f, MOI.VectorConstantChange([2, 2]))
            @test MOI.constant(f) == [2, 2]
            g = deepcopy(f)
            @test g ≈ f
            f = MOIU.modify_function(f, MOI.MultirowChange(z, [(1, 3)]))
            f = MOIU.modify_function(f, StochasticPrograms.DecisionMultirowChange(y, [(2, 3)]))
            @test !(g ≈ f)
            @test g.variable_part.terms == MOI.VectorAffineTerm{Int}[]
            @test g.decision_part.terms == MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([2, 4], [x, y]))
            @test f.variable_part.terms == [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(3, z))]
            @test f.decision_part.terms == MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([2, 3], [x, y]))
            f = MOIU.modify_function(f, StochasticPrograms.DecisionMultirowChange(x, [(1, 0), (2, 1)]))
            @test f.decision_part.terms == MOI.VectorAffineTerm.([2, 2], MOI.ScalarAffineTerm.([3, 1], [y, x]))
        end
    end
end

function runtests()
    x = MOI.VariableIndex(1)
    fx = SingleDecision(x)
    y = MOI.VariableIndex(2)
    fy = SingleDecision(y)
    z = MOI.VariableIndex(3)
    fz = MOI.SingleVariable(z)
    w = MOI.VariableIndex(4)
    fw = MOI.SingleVariable(w)

    @testset "DecisionFunctions" begin
        for name in names(@__MODULE__; all = true)
            if !startswith("$(name)", "test_")
                continue
            end
            f = getfield(@__MODULE__, name)
            @testset "$(name)" begin
                f(x, fx, y, fy, z, fz, w, fw)
            end
        end
    end
end

end


text:
# frexp,ldexp,significand,exponent
@test frexp(12.8) == (0.8,4)
@test ldexp(0.8,4) == 12.8
@test significand(12.8) == 1.6
@test exponent(12.8) == 3

# degree-based trig functions
for T = (Float32,Float64)
    for x = -400:40:400
        @test_approx_eq_eps sind(convert(T,x))::T convert(T,sin(pi/180*x)) eps(deg2rad(convert(T,x)))
        @test_approx_eq_eps cosd(convert(T,x))::T convert(T,cos(pi/180*x)) eps(deg2rad(convert(T,x)))
    end
    for x = 0.0:180:720
        @test sind(convert(T,x)) === zero(T)
        @test sind(-convert(T,x)) === -zero(T)
    end

    for x = -3:0.3:3
        @test_approx_eq_eps sinpi(convert(T,x))::T convert(T,sin(pi*x)) eps(pi*convert(T,x))
        @test_approx_eq_eps cospi(convert(T,x))::T convert(T,cos(pi*x)) eps(pi*convert(T,x))
    end
    for x = 0.0:1.0:4.0
        @test sinpi(convert(T,x)) === zero(T)
        @test sinpi(-convert(T,x)) === -zero(T)
    end
end

# check type stability
for T = (Float32,Float64,BigFloat)
    for f = (sind,cosd,sinpi,cospi)
        @test Base.return_types(f,(T,)) == [T]
    end
end

# error functions
@test_approx_eq erf(1) 0.84270079294971486934
@test_approx_eq erfc(1) 0.15729920705028513066
@test_approx_eq erfcx(1) 0.42758357615580700442
@test_approx_eq erfi(1) 1.6504257587975428760
@test_approx_eq erfinv(0.84270079294971486934) 1
@test_approx_eq erfcinv(0.15729920705028513066) 1
@test_approx_eq dawson(1) 0.53807950691276841914
# TODO: complex versions only supported on 64-bit for now
@unix_only if WORD_SIZE==64
    @test_approx_eq erf(1+2im) -0.53664356577856503399-5.0491437034470346695im
    @test_approx_eq erfc(1+2im) 1.5366435657785650340+5.0491437034470346695im
    @test_approx_eq erfcx(1+2im) 0.14023958136627794370-0.22221344017989910261im
    @test_approx_eq erfi(1+2im) -0.011259006028815025076+1.0036063427256517509im
    @test_approx_eq dawson(1+2im) -13.388927316482919244-11.828715103889593303im
end
for x in logspace(-200, -0.01)
    @test_approx_eq_eps erf(erfinv(x)) x 1e-12*x
    @test_approx_eq_eps erf(erfinv(-x)) -x 1e-12*x
    @test_approx_eq_eps erfc(erfcinv(2*x)) 2*x 1e-12*x
    if x > 1e-20
        xf = float32(x)
        @test_approx_eq_eps erf(erfinv(xf)) xf 1e-5*xf
        @test_approx_eq_eps erf(erfinv(-xf)) -xf 1e-5*xf
        @test_approx_eq_eps erfc(erfcinv(2xf)) 2xf 1e-5*xf
    end
end

# airy
@test_approx_eq airy(1.8) 0.0470362168668458052247
@test_approx_eq airyprime(1.8) -0.0685247801186109345638
@test_approx_eq airybi(1.8) 2.595869356743906290060
@test_approx_eq airybiprime(1.8) 2.98554005084659907283
@test_throws Base.Math.AmosException airy(200im)
@test_throws Base.Math.AmosException airybi(200)

# besselh
true_h133 = 0.30906272225525164362 - 0.53854161610503161800im
@test_approx_eq besselh(3,1,3) true_h133
@test_approx_eq besselh(-3,1,3) -true_h133
@test_approx_eq besselh(3,2,3) conj(true_h133)
@test_approx_eq besselh(-3,2,3) -conj(true_h133)
@test_throws Base.Math.AmosException besselh(1,0)

# besseli
true_i33 = 0.95975362949600785698
@test_approx_eq besseli(3,3) true_i33
@test_approx_eq besseli(-3,3) true_i33
@test_approx_eq besseli(3,-3) -true_i33
@test_approx_eq besseli(-3,-3) -true_i33
@test_throws Base.Math.AmosException besseli(1,1000)

# besselj
@test besselj(0,0) == 1
for i = 1:5
    @test besselj(i,0) == 0
    @test besselj(-i,0) == 0
end

j33 = besselj(3,3.)
@test besselj(3,3) == j33
@test besselj(-3,-3) == j33
@test besselj(-3,3) == -j33
@test besselj(3,-3) == -j33

j43 = besselj(4,3.)
@test besselj(4,3) == j43
@test besselj(-4,-3) == j43
@test besselj(-4,3) == j43
@test besselj(4,-3) == j43

@test_approx_eq j33 0.30906272225525164362
@test_approx_eq j43 0.13203418392461221033
@test_throws DomainError    besselj(0.1, -0.4)
@test_approx_eq besselj(0.1, complex(-0.4)) 0.820421842809028916 + 0.266571215948350899im
@test_approx_eq besselj(3.2, 1.3+0.6im) 0.01135309305831220201 + 0.03927719044393515275im
@test_approx_eq besselj(1, 3im) 3.953370217402609396im
@test_throws Base.Math.AmosException besselj(20,1000im)

# besselk
true_k33 = 0.12217037575718356792
@test_approx_eq besselk(3,3) true_k33
@test_approx_eq besselk(-3,3) true_k33
true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im
@test_throws DomainError besselk(3,-3)
@test_approx_eq besselk(3,complex(-3)) true_k3m3
@test_approx_eq besselk(-3,complex(-3)) true_k3m3
@test_throws Base.Math.AmosException besselk(200,0.01)
# issue #6564
@test besselk(1.0,0.0) == Inf

# bessely
y33 = bessely(3,3.)
@test bessely(3,3) == y33
@test_approx_eq bessely(-3,3) -y33
@test_approx_eq y33 -0.53854161610503161800
@test_throws DomainError bessely(3,-3)
@test_approx_eq bessely(3,complex(-3)) 0.53854161610503161800 - 0.61812544451050328724im
@test_throws Base.Math.AmosException bessely(200.5,0.1)

# issue #6653
for f in (besselj,bessely,besseli,besselk,hankelh1,hankelh2)
    @test_approx_eq f(0,1) f(0,complex128(1))
    @test_approx_eq f(0,1) f(0,complex64(1))
end

# beta, lbeta
@test_approx_eq beta(3/2,7/2) 5π/128
@test_approx_eq beta(3,5) 1/105
@test_approx_eq lbeta(5,4) log(beta(5,4))
@test_approx_eq beta(5,4) beta(4,5)
@test_approx_eq beta(-1/2, 3) -16/3
@test_approx_eq lbeta(-1/2, 3) log(16/3)

# gamma, lgamma (complex argument)
@test gamma(1:25) == gamma(Float64[1:25])
@test_approx_eq gamma(1/2) sqrt(π)
@test_approx_eq gamma(-1/2) -2sqrt(π)
@test_approx_eq lgamma(-1/2) log(abs(gamma(-1/2)))
@test_approx_eq lgamma(1.4+3.7im) -3.7094025330996841898 + 2.4568090502768651184im
@test_approx_eq lgamma(1.4+3.7im) log(gamma(1.4+3.7im))

# digamma
for elty in (Float32, Float64)

    @test_approx_eq digamma(convert(elty, 9)) convert(elty, 2.140641477955609996536345)
    @test_approx_eq digamma(convert(elty, 2.5)) convert(elty, 0.7031566406452431872257)
    @test_approx_eq digamma(convert(elty, 0.1)) convert(elty, -10.42375494041107679516822)
    @test_approx_eq digamma(convert(elty, 7e-4)) convert(elty, -1429.147493371120205005198)
    @test_approx_eq digamma(convert(elty, 7e-5)) convert(elty, -14286.29138623969227538398)
    @test_approx_eq digamma(convert(elty, 7e-6)) convert(elty, -142857.7200612932791081972)
    @test_approx_eq digamma(convert(elty, 2e-6)) convert(elty, -500000.5772123750382073831)
    @test_approx_eq digamma(convert(elty, 1e-6)) convert(elty, -1000000.577214019968668068)
    @test_approx_eq digamma(convert(elty, 7e-7)) convert(elty, -1428572.005785942019703646)
    @test_approx_eq digamma(convert(elty, -0.5)) convert(elty, .03648997397857652055902367)
    @test_approx_eq digamma(convert(elty, -1.1)) convert(elty,  10.15416395914385769902271)

    @test_approx_eq digamma(convert(elty, 0.1)) convert(elty, -10.42375494041108)
    @test_approx_eq digamma(convert(elty, 1/2)) convert(elty, -γ - log(4))
    @test_approx_eq digamma(convert(elty, 1)) convert(elty, -γ)
    @test_approx_eq digamma(convert(elty, 2)) convert(elty, 1 - γ)
    @test_approx_eq digamma(convert(elty, 3)) convert(elty, 3/2 - γ)
    @test_approx_eq digamma(convert(elty, 4)) convert(elty, 11/6 - γ)
    @test_approx_eq digamma(convert(elty, 5)) convert(elty, 25/12 - γ)
    @test_approx_eq digamma(convert(elty, 10)) convert(elty, 7129/2520 - γ)
end

# trigamma
for elty in (Float32, Float64)
    @test_approx_eq trigamma(convert(elty, 0.1)) convert(elty, 101.433299150792758817)
    @test_approx_eq trigamma(convert(elty, 1/2)) convert(elty, π^2/2)
    @test_approx_eq trigamma(convert(elty, 1)) convert(elty, π^2/6)
    @test_approx_eq trigamma(convert(elty, 2)) convert(elty, π^2/6 - 1)
    @test_approx_eq trigamma(convert(elty, 3)) convert(elty, π^2/6 - 5/4)
    @test_approx_eq trigamma(convert(elty, 4)) convert(elty, π^2/6 - 49/36)
    @test_approx_eq trigamma(convert(elty, 5)) convert(elty, π^2/6 - 205/144)
    @test_approx_eq trigamma(convert(elty, 10)) convert(elty, π^2/6 - 9778141/6350400)
end

# invdigamma
for elty in (Float32, Float64)
    for val in [0.001, 0.01, 0.1, 1.0, 10.0]
        @test abs(invdigamma(digamma(convert(elty, val))) - convert(elty, val)) < 1e-8
    end
end

@test_approx_eq polygamma(20, 7.) -4.644616027240543262561198814998587152547

# eta, zeta
@test_approx_eq eta(1) log(2)
@test_approx_eq eta(2) pi^2/12
@test_approx_eq zeta(0) -0.5
@test_approx_eq zeta(2) pi^2/6
@test_approx_eq zeta(4) pi^4/90

# quadgk
@test_approx_eq quadgk(cos, 0,0.7,1)[1] sin(1)
@test_approx_eq quadgk(x -> exp(im*x), 0,0.7,1)[1] (exp(1im)-1)/im
@test_approx_eq quadgk(x -> exp(im*x), 0,1im)[1] -1im*expm1(-1)
@test_approx_eq_eps quadgk(cos, 0,BigFloat(1),order=40)[1] sin(BigFloat(1)) 1000*eps(BigFloat)
@test_approx_eq quadgk(x -> exp(-x), 0,0.7,Inf)[1] 1.0
@test_approx_eq quadgk(x -> exp(x), -Inf,0)[1] 1.0
@test_approx_eq quadgk(x -> exp(-x^2), -Inf,Inf)[1] sqrt(pi)
@test_approx_eq quadgk(x -> [exp(-x), exp(-2x)], 0, Inf)[1] [1,0.5]
@test_approx_eq quadgk(cos, 0,0.7,1, norm=abs)[1] sin(1)

# Ensure subnormal flags functions don't segfault
@test any(ccall("jl_zero_subnormals", Uint8, (Uint8,), 1) .== [0x00 0x01])
@test any(ccall("jl_zero_subnormals", Uint8, (Uint8,), 0) .== [0x00 0x01])

# isqrt (issue #4884)
@test isqrt(9223372030926249000) == 3037000498
@test isqrt(typemax(Int128)) == int128("13043817825332782212")
@test isqrt(int128(typemax(Int64))^2-1) == 9223372036854775806
@test isqrt(0) == 0
for i = 1:1000
    n = rand(Uint128)
    s = isqrt(n)
    @test s*s <= n
    @test (s+1)*(s+1) > n
    n = rand(Uint64)
    s = isqrt(n)
    @test s*s <= n
    @test (s+1)*(s+1) > n
end

ext: jl

max_stars_repo_name: vchuravy/julia

max_issues_repo_name: vchuravy/julia

max_forks_repo_name: vchuravy/julia


text:
# Types for Boundary Conditions
# Abstracts
abstract type BoundaryCondition end
abstract type HomogeneousBoundaryCondition <: BoundaryCondition end
abstract type NonhomogeneousBoundaryCondition <: BoundaryCondition end

# Types for Differences
abstract type DiscretizationMethod end
struct DefaultFirstDifference <: DiscretizationMethod end
struct ForwardFirstDifference <: DiscretizationMethod end
struct BackwardFirstDifference <: DiscretizationMethod end
struct CentralSecondDifference <: DiscretizationMethod end
struct JumpProcess{T} <: DiscretizationMethod where T <: AbstractArray
    jumps::T
end

"""
JumpProcess(x̄, jumps::AbstractArray, truncate = (:interior, :interior))

Returns a DiscretizationMethod object that can be used to construct
a discretized operator jump process.

`jumps` is a `(length(x̄)-2)`-vector whose ith element is an integer that represents
a jump size in index and direction (by sign) from `i`th element of `interiornodes(x̄)` and
the first and second elements of `truncate` represent truncation location for
the lower bound and upper bound when the jump is out of the truncated boundary
(`:interior` for the first/last element of `interiornodes(̄x)` and `:exterior` for
the first/last element of `x̄`). The default parameter is `(:interior, :interior)`.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> JumpProcess(x̄, [-1; -1; -1; -1], (:interior, :interior))
JumpProcess{Array{Int64,1}}([0, -1, -1, -1])
```
"""
function JumpProcess(x̄::AbstractArray, jumps::AbstractArray,
                    truncate = (:interior, :interior))
    # each ith element of jumps defines the jump direction (negative/positive) and size
    # where i ranges from 1 to M (interior points)
    @assert length(x̄) == length(jumps) + 2
    M = length(x̄) - 2

    # define default destinations and truncation indices
    destinations = Array(1:M) + jumps
    cutoff_lb = (truncate[1] == :interior) ? 1 : 0
    cutoff_ub = (truncate[2] == :interior) ? M : M+1

    # perform truncations
    destinations[destinations .< cutoff_lb] .= cutoff_lb
    destinations[destinations .> cutoff_ub] .= cutoff_ub
    jumps_after_truncation = Array{Int64}(destinations - Array(1:M))

    return JumpProcess(jumps_after_truncation)
end

"""
JumpProcess(x̄, uniform_jump::AbstractArray, truncate = (:interior, :interior))

Returns a DiscretizationMethod object that can be used to construct
a discretized operator jump process.

`uniform_jump` is a scalar Int64 that represents a jump size in index
and direction (by sign) from all elements of `interiornodes(x̄)` and
the first and second elements of `truncate` represent truncation location for
the lower bound and upper bound when the jump is out of the truncated boundary
(`:interior` for the first/last element of `interiornodes(̄x)` and `:exterior` for
the first/last element of `x̄`). The default parameter is `(:interior, :interior)`.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> JumpProcess(x̄, -1)
JumpProcess{Array{Int64,1}}([0, -1, -1, -1])
```
"""
JumpProcess(x̄::AbstractArray, uniform_jump::Int64, truncate = (:interior, :interior)) = JumpProcess(x̄, uniform_jump*ones(length(x̄) - 2), truncate)

"""
JumpProcess(x̄, jumpf::Function, truncate = (:interior, :interior))

Returns a DiscretizationMethod object that can be used to construct
a discretized operator jump process.

`jumpf` is a function that takes an element of `interiornodes(x̄)` and returns
a scalar that represents the corresponding nominal jump size and direction (by sign) and
the first and second elements of `truncate` represent truncation location for
the lower bound and upper bound when the jump is out of the truncated boundary
(`:interior` for the first/last element of `interiornodes(̄x)` and `:exterior` for
the first/last element of `x̄`). The default parameter is `(:interior, :interior)`.
The code uses nearest-neighbour rule to determine the indices of destinations according
to `jumpf.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> jumpf(x) = -1.4
jumpf (generic function with 1 method)

julia> JumpProcess(x̄, jumpf)
JumpProcess{Array{Int64,1}}([0, -1, -1, -1])
```
"""
function JumpProcess(x̄::AbstractArray, jumpf::Function,
    truncate = (:interior, :interior))
    # each ith element of jumps defines the jump direction (negative/positive) and size
    # where i ranges from 1 to M (interior points)
    x = interiornodes(x̄)
    M = length(x)

    # define default destinations and truncation indices
    destinations = x + jumpf.(x)

    cutoff_lb = (truncate[1] == :interior) ? x[1] : x̄[1]
    cutoff_ub = (truncate[2] == :interior) ? x[end] : x̄[end]

    # perform truncations
    destinations[destinations .< cutoff_lb] .= cutoff_lb
    destinations[destinations .> cutoff_ub] .= cutoff_ub

    # find the corresponding destination indices and jump indices
    destinations_indices_on_x̄ = (destination -> findnearestindex(x̄, destination)).(destinations)
    current_indices_on_x̄ = Array(2:(M+1))
    jumps_after_truncation = destinations_indices_on_x̄ - current_indices_on_x̄

    return JumpProcess(jumps_after_truncation)
end

"""
JumpProcess(x̄, uniform_jump_size::Real, truncate = (:interior, :interior))

Returns a DiscretizationMethod object that can be used to construct
a discretized operator jump process.

`uniform_jump_size` is a scalar in Real that represents a nominal jump size
and direction (by sign) from all elements of `interiornodes(x̄)` and
the first and second elements of `truncate` represent truncation location for
the lower bound and upper bound when the jump is out of the truncated boundary
(`:interior` for the first/last element of `interiornodes(̄x)` and `:exterior` for
the first/last element of `x̄`). The default parameter is `(:interior, :interior)`.
The code uses nearest-neighbour rule to determine the indices of destinations according
to `jumpf.
# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> JumpProcess(x̄, -1.4)
JumpProcess{Array{Int64,1}}([0, -1, -1, -1])
```
"""
JumpProcess(x̄::AbstractArray, uniform_jump_size::Real, truncate = (:interior, :interior)) = JumpProcess(x̄, (x -> uniform_jump_size), truncate)

# Concretes
struct Reflecting <: HomogeneousBoundaryCondition end
@with_kw struct Mixed{T} <: HomogeneousBoundaryCondition where T <: Real
    ξ::T = 0.0
    direction = :auto # :forward/:backward/:auto
end
@with_kw struct Absorbing{T} <: HomogeneousBoundaryCondition where T <: Int64
    loc::T = 0 # if loc > 0, apply f(x̄ₘ) = 0 for m <= loc, where x̄₀ is the lower boundary
end
@with_kw struct NonhomogeneousAbsorbing{T} <: NonhomogeneousBoundaryCondition where T <: Real
    S::T  = 0.0
end


text:
"""
    function addUpdateLocal_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions related with Global Updates.
"""
function addUpdateGlobal_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    if "UpdateGlobal" in keys(abm.declaredUpdates)

        #Check updated
        up = symbols_(abm,abm.declaredUpdates["UpdateGlobal"])
        up = up[Bool.((up[:,"placeDeclaration"].==:Model) .* Bool.((up[:,"assigned"].==true) .+ (up[:,"updated"].==true))),:]

        #Construct functions
        f = simpleFirstLoopWrapInFunction_(platform,:globStep_!,
                        :(begin
                            if ic1_ == 1
                                $(abm.declaredUpdates["UpdateGlobal"])
                            end
                        end)
                        )
        f = vectorize_(abm,f,p)

        push!(p.declareF.args,
            f)

        push!(p.execInloop.args,
                :(@platformAdapt globStep_!(ARGS_))
            )
    end

    return nothing
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
function interpolate1D(f::AbstractVector, x::AbstractVector, x0)
    """Returns value of function f at x0 by linear interpolation of provided datapoints.
    x should be ascending. Interpolation is periodic"""

    bounds = getBoundsIndices(x, x0)
    i = bounds[1]
    j = bounds[2]

    # Interpolate and return
    f[i] + (f[j] - f[i])*(x0 - x[i])/(x[j] - x[i])
end

function interpolate2D(f::Matrix, x::AbstractVector, y::AbstractVector, x0, y0)
    """
    Returns value of function f at (x0,y0) by repeated linear interpolation of provided datapoints.
    x should ascend as row index increases. y should ascend as column index increases.
    Interpolation is periodic.
    """

    # Find upper and lower y indices of grid square containing (x0, y0)
    # We don't need the x indices
    bounds = getBoundsIndices(y, y0)
    i = bounds[1]
    j = bounds[2]

    # Intepolate for x at borders of grid (fixed y)
    fi = interpolate1D(f[:,i], x, x0)
    fj = interpolate1D(f[:,j], x, x0)

    return interpolate1D([fi, fj], [y[i], y[j]], y0)
end

function interpolate3D(f::Array{T, 3}, x::AbstractVector, y::AbstractVector, z::AbstractVector, x0, y0, z0) where T
    """Returns value of function f at (x0,y0,z0) by bilinear interpolation of provided datapoints.
    x should be ascending. Interpolation is periodic."""

    # Find bounding incides of grid cube containing (x0, y0, z0)\
    # We don't need the x indices
    i,j = getBoundsIndices(z, z0)

    # Interpolate for f(x0,y0) holding z fixed
    fi = interpolate2D(f[:,:,i], x, y, x0, y0)
    fj = interpolate2D(f[:,:,j], x, y, x0, y0)

    # Interpolate between the planes
    return interpolate1D([fi, fj], [z[i], z[j]], z0)
end

function getBoundsIndices(x::AbstractVector, x0)::Tuple
    """Returns indices i and j of points such that x[i] < x <= x[j].
    Values of x0 > max are treated periodically. Assumes x is ascending."""

    # Make sure x0 is in range
    if minimum(x) > x0 || x0 >= maximum(x)
        @debug "Point out of range of input vectors"
        x0 = rem(x0, maximum(x))
    end

    # Get upper bounds
    j = findfirst(x .>= x0)
    i = j == 1 ? length(x) : j-1

    return (i, j)
end

# This all works, but there's no point in loading it in production
# using Plots
# function interpolate2DTest()
#     data = rand(Float64, (4,4))
#     x = Vector(1:4)
#     y = Vector(1:4)

#     xpoints = 1:0.1:4
#     ypoints = 1:0.1:4

#     int2d = [interpolate2D(data, x, y, x0, y0) for x0 in xpoints, y0 in ypoints]
#     display(surface(xpoints, ypoints, int2d))
# end

# function interpolate3DTest()

#     data = rand(Float64, (4,4,4))
#     x = Vector(1:4)
#     y = Vector(1:4)
#     z = Vector(1:4)

#     xpoints = 1:0.1:4
#     ypoints = 1:0.1:4
#     zpoints = 1:0.1:4

#     [interpolate3D(data, x, y, z, x0, y0, z0) for x0 in xpoints, y0 in ypoints, z0 in zpoints]
# end
# interpolate2DTest()
# interpolate3DTest()


text:
using StaticArrays
using StaticArrays.ImmutableArrays
using PyPlot
export Polygon, isInside, isIntersect

const Vec2f = SVector{2, Float64}
const Vec4f = SVector{4, Float64}

abstract type Polygonic end

struct Polygon <: Polygonic
    N::Int
    V::Vector{Vec2f}
    function Polygon(V_::AbstractVector) #V nust be counterclockwise
        N = length(V_)
        V = Vec2f[]
        for v in V_
            push!(V, Vec2f(v))
        end
        new(N, V)
    end
end

struct Rectangle <: Polygonic
    N::Int
    V::Vector{Vec2f}
    function Rectangle(center::T, width::Float64, height::Float64) where {T<:AbstractVector}
        V1 = Vec2f(center + [-width, -height]*0.5)
        V2 = Vec2f(center + [+width, -height]*0.5)
        V3 = Vec2f(center + [+width, +height]*0.5)
        V4 = Vec2f(center + [-width, +height]*0.5)
        new(4, [V1, V2, V3, V4])
    end
end

@inline function isInside(this::T, q::Vec2f) where {T<:Polygonic}
    # note: n'x = n'b
    @inbounds q_vec = [q[1]; q[2]]
    for n = 1:this.N
        @inbounds p1 = this.V[n]
        @inbounds n<this.N ? p2 = this.V[n+1] : p2 = this.V[1]
        @inbounds u = [p2[1]-p1[1]; p2[2]-p1[2]] #p2-p1
        @inbounds n = [0 1; -1 0]*u # normla vector
        @inbounds px = [q[1]-p1[1]; q[2]-p1[2]] # x-p1
        px'*n>0 && return false
    end
    return true
end

@inline function isIntersect(this::T, q1_::AbstractVector, q2_::AbstractVector) where {T<:Polygonic}
    @views q1 = Vec2f(q1_[1:2])
    @views q2 = Vec2f(q2_[1:2])
    return isIntersect(this, q1, q2)
end

@inline function isIntersect(this::T, q1::Vec2f, q2::Vec2f)::Bool where {T<:Polygonic}
    for n = 1:this.N
        p1 = this.V[n]
        if n<this.N
            p2=this.V[n+1]
        else
            p2=this.V[1]
        end
        isIntersect(p1, p2, q1, q2) && return true
    end
    return false
end

@inline function isIntersect(this::T, q_seq::Vector{Vec2f})::Bool where {T<:Polygonic}
    for q in q_seq
        isInside(this, q) && return true
    end
    return false
end

function show(this::T, color=:red) where {T<:Polygonic}
    for n = 1:this.N
        if n<this.N
            x = [this.V[n][1], this.V[n+1][1]]
            y = [this.V[n][2], this.V[n+1][2]]
        else
            x = [this.V[n][1], this.V[1][1]]
            y = [this.V[n][2], this.V[1][2]]
        end
        plot(x, y, "r-")
    end
end

@inline function isIntersect(p1::Vec2f, p2::Vec2f, q1::Vec2f, q2::Vec2f)::Bool
    # solve [p2-p1, q2-q1]*[s; t]=[q1-p1] or A*[s; t]=B
    # for computational efficiency, not using inv() func
    # A = [a, b; c, d]
    # B = [e; f]
    a = p2[1]-p1[1]
    b = -(q2[1]-q1[1])
    c = p2[2]-p1[2]
    d = -(q2[2]-q1[2])
    e = q1[1]-p1[1]
    f = q1[2]-p1[2]
    det = a*d-b*c
    abs(det)<1e-5 && return false

    s = 1/det*(d*e-b*f)
    t = 1/det*(-c*e+a*f)
    return (0.0<=s<=1.0) && (0.0<=t<=1.0)
end

function test()
    """
    v1 = [0.0, 0.0]
    v2 = [1.0, 0.0]
    v3 = [1.0, 1.0]
    v4 = [0.0, 1.0]
    V = [v1, v2, v3, v4]
    P = Polygon(V)
    """
    P = Rectangle([0.5 0.5], 1.0, 1.0)
    isInside(P, Vec2f([1.2, 0.8]))
end


text:
@testset "runtestCountdown" begin
    println("Testing runtestCountdown")
    include("../src/bfgaCountdown.jl")
    using .bfgaCountdown

    include("../src/Bf.jl")
    using .BfInterpreter

    @time model = test_serial()
    prems = model.refSet[length(model.refSet)]
    println("Generation : $(model.params.currentGeneration)")
    println("& fitness : $(fitness(prems, model.instructionsSet))")
    #bfcode :: Array{Float64,1}  = prems.dna
    println("code : $(prems.program) " )# $(join(genesToBfInstructions(bfcode), "") ) ")

    #_trainingExamples = [ "00" , "01" , "10" , "11" ]
    #_trainingResults = [ '1', '0', '0', '1' ]

    #bft = bfType(bfcode)
    for t in _trainingExamples
        output, _ = execute(prems.program, t, model.instructionsSet)
        @test length(output) >= (t[1] +1)
        for i in 0:t[1]
            @test t[1]-i == Int(output[(i+1)])
        end
    end

    # do the test_parallel after ...
    #test_parallel()
end


text:
using LinearAlgebra

function cheb(N)
    if N==0
        D=0; x=1
        return D,x
    else
        x= cos.(pi * collect(0:N) ./ N)
        c=[2; ones(N-1, 1); 2] .* (-1).^collect(0:N);
        X = repeat(x, 1, N+1);
        dX=X .- X';
        D = (c .* (1 ./ c)')  ./ (dX .+ Diagonal( ones(N+1) ) );
        D = D .- Diagonal( vec( sum(D', dims=1) ) );
        return D,x
    end
end

function clencurt(N)

     θ = π * collect(0:N) ./N
     x = cos.(θ)
     w = zeros(1, N+1)
    ii = collect(2:N)
     v = ones(N-1,1)

    if mod(N,2) == 0
          w[1] = 1/(N^2 - 1)
        w[N+1] = w[1]
        for k=1:N/2-1
             v = v .- 2 * cos.(2 * k * θ[ii])/(4 * k^2 - 1)
        end
        v = v - cos.(N*θ[ii]) / (N^2 - 1)
    else
          w[1] = 1/N^2
        w[N+1] = w[1]
        for k=1:(N-1)/2
            v = v .- 2 * cos.(2 * k * θ[ii])/(4 * k^2 - 1)
        end
    end

    w[ii] = 2*v/N

    return x, w
end


text:
struct FreeFlow <: AbstractFlowDynamics end

function collision_time(b::MovableBall, Π::FixedPlane, ::FreeFlow)
    @unpack r, x, v = b
    @unpack n, p = Π

    if v ⋅ n < 0  # travelling in wrong direction
        return Inf
    end

    return (-r - (x - p)⋅n) / (v ⋅ n)
end

function collision_time(b1::MovableBall, b2::MovableBall, ::FreeFlow)
    Δx = b1.x - b2.x
    Δv = b1.v - b2.v

    b = Δx⋅Δv

    b > 0 && return Inf  # moving away from one another

    a = normsq(Δv)
    c = normsq(Δx) - (b1.r + b2.r)^2

    discriminant = b^2 - a*c

    if discriminant ≥ 0
        d = √discriminant

        # t1 = (-b + d) / a
        t2 = (-b - d) / a

        if t2 > 0
            return t2
        end
    end

    return Inf   # no collision
end

#
# function find_collision(b::MovableBall{N,T}, box::RectangularBox{N,T}, ::FreeFlow) where {N,T}
#
#     walls = box.walls
#
#     tol = 1e-10
#
#     min_collision_time = Inf
#     which = -1
#
#     for i in 1:length(walls)
#         t = collision_time(b, walls[i])
#         if tol < t < min_collision_time
#             min_collision_time = t
#             which = i
#         end
#     end
#
#     return (which, min_collision_time)
# end

"Flow ball for a time t"
function flow!(b::MovableBall, t, ::FreeFlow)
    b.x += b.v * t
end

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
@testset "MRCHeader" begin
    @testset "MRCHeader()" begin
        h = MRCHeader()
        @test h.nx === Int32(0)
        @test h.ny === Int32(0)
        @test h.nz === Int32(0)
        @test h.mode === Int32(0)
        @test h.nxstart === Int32(0)
        @test h.nystart === Int32(0)
        @test h.nzstart === Int32(0)
        @test h.mx === Int32(0)
        @test h.my === Int32(0)
        @test h.mz === Int32(0)
        @test h.cella_x === Float32(0)
        @test h.cella_y === Float32(0)
        @test h.cella_z === Float32(0)
        @test h.cellb_alpha === Float32(90)
        @test h.cellb_beta === Float32(90)
        @test h.cellb_gamma === Float32(90)
        @test h.mapc === Int32(1)
        @test h.mapr === Int32(2)
        @test h.maps === Int32(3)
        @test h.dmin > h.dmax # test undefined
        @test h.dmean > h.dmax || h.dmean < h.dmin # test undefined
        @test h.ispg === Int32(1)
        @test h.nsymbt === Int32(0)
        @test h.extra1 === Tuple(zeros(UInt8, 8))
        @test h.exttyp === ""
        @test h.nversion === Int32(20140)
        @test h.extra2 === Tuple(zeros(UInt8, 84))
        @test h.origin_x === Float32(0)
        @test h.origin_y === Float32(0)
        @test h.origin_z === Float32(0)
        @test h.map === "MAP "
        @test h.machst === Tuple(MRC.machstfrombyteorder())
        @test h.rms === Float32(-1)
        @test h.nlabl === Int32(1)
        @test length(h.label) == 10
        @test length(h.label[1]) < 80
        @test startswith(h.label[1], "Created by MRC.jl at ")
    end

    @testset "MRCHeader(kwargs...)" begin
        h = MRCHeader(; nx=10)
        @test h.nx === Int32(10)
        h2 = MRCHeader(; nsymbt=3)
        @test h2.nsymbt === Int32(3)
        h3 = MRCHeader(; ispg=2)
        @test h3.ispg === Int32(2)
    end
end

@testset "show(::MRCHeader)" begin
    h = MRCHeader()
    @test startswith(
        sprint(show, "text/plain", h),
        """
MRCHeader(
    nx = 0,
    ny = 0,
    nz = 0,
    mode = 0,
    nxstart = 0,
    nystart = 0,
    nzstart = 0,
    mx = 0,
    my = 0,
    mz = 0,
    cella_x = 0.0f0,
    cella_y = 0.0f0,
    cella_z = 0.0f0,
    cellb_alpha = 90.0f0,
    cellb_beta = 90.0f0,
    cellb_gamma = 90.0f0,
    mapc = 1,
    mapr = 2,
    maps = 3,
    dmin = 0.0f0,
    dmax = -1.0f0,
    dmean = -2.0f0,
    ispg = 1,
    nsymbt = 0,
    extra1 = (0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
    exttyp = "",
    nversion = 20140,
    extra2 = (0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
    origin_x = 0.0f0,
    origin_y = 0.0f0,
    origin_z = 0.0f0,
    map = "MAP ",
    machst = (0x44, 0x44, 0x00, 0x00),
    rms = -1.0f0,
    nlabl = 1,
""",
    )
end

@testset "setproperty!(::MRCHeader)" begin
    h = MRCHeader()
    @test h.nx == 0
    h.nx = 10
    @test h.nx == 10
end

@testset "copyto!(::MRCHeader, ::MRCHeader)" begin
    h = MRCHeader()
    h2 = MRCHeader(; nx=10, ny=10, nz=10)
    copyto!(h, h2)
    @test h.nx == 10
    @test h.ny == 10
    @test h.nz == 10
end

@testset "sizeoffield($name, $type)" for (name, type, s) in (
    (:exttyp, String, 4),
    (:map, String, 4),
    (:label, Vector{String}, 800),
    (:testfloat, Float32, sizeof(Float32)),
    (:testint, Int32, sizeof(Int32)),
    (:testuintbool, NTuple{8,UInt8}, sizeof(NTuple{8,UInt8})),
)
    @test MRC.sizeoffield(name, type) == s
end

@testset "entrytobytes" begin
    @test MRC.entrytobytes(:exttyp, "abc") == [0x61, 0x62, 0x63, 0x00]
    @test MRC.entrytobytes(:exttyp, "abcde") == [0x61, 0x62, 0x63, 0x64]
    @test MRC.entrytobytes(:map, "abc") == [0x61, 0x62, 0x63, 0x00]
    @test MRC.entrytobytes(:map, "abcde") == [0x61, 0x62, 0x63, 0x64]
    @test MRC.entrytobytes(:label, ("abcd", ntuple(_ -> "", 9)...))[1:80] ==
          [0x61; 0x62; 0x63; 0x64; zeros(UInt8, 76)]
    @test MRC.entrytobytes(:testfloat, Float32(3)) == reinterpret(UInt8, [Float32(3)])
    @test MRC.entrytobytes(:testint, Int32(3)) == reinterpret(UInt8, [Int32(3)])
    @test MRC.entrytobytes(:testuintbool, (0x01, 0x02, 0x03, 0x04)) ==
          [0x01, 0x02, 0x03, 0x04]
end

@testset "fieldoffsets" begin
    @test MRC.fieldoffsets(MRCHeader) == [
        0,
        4,
        8,
        12,
        16,
        20,
        24,
        28,
        32,
        36,
        40,
        44,
        48,
        52,
        56,
        60,
        64,
        68,
        72,
        76,
        80,
        84,
        88,
        92,
        96,
        104,
        108,
        112,
        196,
        200,
        204,
        208,
        212,
        216,
        220,
        224,
    ]
end

@testset "size(::MRCHeader)" begin
    h = MRCHeader(; nx=10, ny=20, nz=30)
    @test size(h) == (10, 20, 30)
end

@testset "length(::MRCHeader)" begin
    h = MRCHeader(; nx=10, ny=20, nz=30)
    @test length(h) == 6000
end

# @test "ndims(::MRCHeader)" begin
#     h = MRCHeader()
# end

@testset "cellangles(::MRCHeader)" begin
    h = MRCHeader(; cellb_alpha=10, cellb_beta=20, cellb_gamma=30)
    @test cellangles(h) == (10.0f0, 20.0f0, 30.0f0)
end

@testset "cellangles!(::MRCHeader)" begin
    h = MRCHeader()
    @test cellangles(h) == (90.0f0, 90.0f0, 90.0f0)
    cellangles!(h, (10, 20, 30))
    @test cellangles(h) == (10.0f0, 20.0f0, 30.0f0)
end

@testset "cellsize(::MRCHeader)" begin
    h = MRCHeader(; cella_x=10, cella_y=20, cella_z=30)
    @test cellsize(h) == (10.0f0, 20.0f0, 30.0f0)
end

@testset "cellsize!(::MRCHeader)" begin
    h = MRCHeader()
    @test cellsize(h) == (0.0f0, 0.0f0, 0.0f0)
    cellsize!(h, (10, 20, 30))
    @test cellsize(h) == (10.0f0, 20.0f0, 30.0f0)
end

@testset "gridsize(::MRCHeader)" begin
    h = MRCHeader(; mx=10, my=20, mz=30)
    @test gridsize(h) == Int32.((10, 20, 30))
end

@testset "gridsize!(::MRCHeader)" begin
    h = MRCHeader()
    @test gridsize(h) == Int32.((0, 0, 0))
    gridsize!(h, (10, 20, 30))
    @test gridsize(h) == Int32.((10, 20, 30))
end

@testset "origin(::MRCHeader)" begin
    h = MRCHeader(; origin_x=10, origin_y=20, origin_z=30)
    @test origin(h) == (10.0f0, 20.0f0, 30.0f0)
end

@testset "origin!(::MRCHeader)" begin
    h = MRCHeader()
    @test origin(h) == (0.0f0, 0.0f0, 0.0f0)
    origin!(h, (10, 20, 30))
    @test origin(h) == (10.0f0, 20.0f0, 30.0f0)
end

@testset "start(::MRCHeader)" begin
    h = MRCHeader(; nxstart=10, nystart=20, nzstart=30)
    @test start(h) == Int32.((10, 20, 30))
end

@testset "start!(::MRCHeader)" begin
    h = MRCHeader()
    @test start(h) == Int32.((0, 0, 0))
    start!(h, (10, 20, 30))
    @test start(h) == Int32.((10, 20, 30))
end

@testset "voxelsize" begin
    @testset "voxelsize(::MRCHeader)" begin
        h = MRCHeader(; cella_x=10, cella_y=20, cella_z=30, mx=2, my=5, mz=6)
        @test voxelsize(h) == (5.0f0, 4.0f0, 5.0f0)
    end

    @testset "voxelsize(::MRCHeader, i)" for i in 1:3
        scella = (:cella_x, :cella_y, :cella_z)[i]
        smx = (:mx, :my, :mz)[i]
        h = MRCHeader(; scella => 10, smx => 2)
        @test voxelsize(h, i) == 5.0f0
    end
end

@testset "voxelsize!" begin
    @testset "voxelsize!(::MRCHeader, s)" begin
        h = MRCHeader(; mx=10, my=20, mz=30)
        voxelsize!(h, (5, 10, 15))
        @test voxelsize(h) == Float32.((5, 10, 15))
        @test cellsize(h) == Float32.((50, 200, 450))
    end

    @testset "voxelsize!(::MRCHeader, s, i)" begin
        h = MRCHeader(; mx=10)
        voxelsize!(h, 5, 1)
        @test voxelsize(h, 1) == 5
        @test h.cella_x == 50
    end
end


text:
include("euler/euler.jl")
using .Primes: prime_numbers

function compute(n::Int)::Int
    primes = prime_numbers(10 ^ n ÷ n)
    consecutive = 0
    i = 10 ^ (n - 1)
    while true
        if i ∈ primes
            consecutive = 0
        else
            divisors = 0
            t = i
            for prime ∈ primes
                if t % prime == 0
                    divisors += 1
                    t ÷= prime
                end
                if prime > t || t == 1 || divisors == n
                    break
                end
            end
            if divisors == n
                consecutive += 1
            else
                consecutive = 0
            end
        end
        if consecutive == n return i - n + 1 end
        i += 1
    end
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler