text:
module TestBridgeOptimizer

using Test

using MathOptInterface

const MOI = MathOptInterface

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

include("utilities.jl")

struct DummyModelAttribute <: MOI.AbstractModelAttribute end
struct DummyEvaluator <: MOI.AbstractNLPEvaluator end
struct DummyVariableAttribute <: MOI.AbstractVariableAttribute end
struct DummyConstraintAttribute <: MOI.AbstractConstraintAttribute end

function test_subsitution_of_variables()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged = MOI.Bridges.Variable.Vectorize{Float64}(mock)
    x, cx = MOI.add_constrained_variable(bridged, MOI.GreaterThan(1.0))
    y = MOI.get(mock, MOI.ListOfVariableIndices())[1]
    c2x = MOI.add_constraint(bridged, 2.0x, MOI.GreaterThan(1.0))
    @test MOI.get(bridged, MOI.ConstraintFunction(), c2x) ≈ 2.0x
    @test MOI.get(bridged, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(1.0)
    # The constant was moved to the set
    @test MOI.get(mock, MOI.ConstraintFunction(), c2x) ≈ 2.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(-1.0)

    MOI.set(bridged, MOI.ConstraintSet(), c2x, MOI.GreaterThan(2.0))
    @test MOI.get(bridged, MOI.ConstraintFunction(), c2x) ≈ 2.0x
    @test MOI.get(bridged, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(2.0)
    @test MOI.get(mock, MOI.ConstraintFunction(), c2x) ≈ 2.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(0.0)

    MOI.set(bridged, MOI.ConstraintFunction(), c2x, 3.0x)
    @test MOI.get(bridged, MOI.ConstraintFunction(), c2x) ≈ 3.0x
    @test MOI.get(bridged, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(2.0)
    @test MOI.get(mock, MOI.ConstraintFunction(), c2x) ≈ 3.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(-1.0)

    MOI.set(bridged, MOI.ConstraintSet(), c2x, MOI.GreaterThan(4.0))
    @test MOI.get(bridged, MOI.ConstraintFunction(), c2x) ≈ 3.0x
    @test MOI.get(bridged, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(4.0)
    @test MOI.get(mock, MOI.ConstraintFunction(), c2x) ≈ 3.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(1.0)

    MOI.set(bridged, MOI.ObjectiveFunction{typeof(2.0x)}(), 2.0x)
    @test MOI.get(mock, MOI.ObjectiveFunction{typeof(2.0y)}()) ≈ 2.0y + 2.0

    MOI.set(bridged, DummyModelAttribute(), 2.0x + 1.0)
    @test MOI.get(bridged, DummyModelAttribute()) ≈ 2.0x + 1.0
    @test MOI.get(mock, DummyModelAttribute()) ≈ 2.0y + 3.0

    z = MOI.add_variable(bridged)
    for (attr, index) in
        [(DummyVariableAttribute(), z), (DummyConstraintAttribute(), c2x)]
        MOI.set(bridged, attr, index, 1.0x)
        @test MOI.get(bridged, attr, index) ≈ 1.0x
        @test MOI.get(mock, attr, index) ≈ 1.0y + 1.0

        MOI.set(bridged, attr, [index], [3.0x + 1.0z])
        @test MOI.get(bridged, attr, [index])[1] ≈ 3.0x + 1.0z
        @test MOI.get(mock, attr, [index])[1] ≈ 3.0y + 1.0z + 3.0
    end
    return
end

function test_raw_solver()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged = MOI.Bridges.Variable.Vectorize{Float64}(mock)
    return MOI.get(bridged, MOI.RawSolver()) === mock
end

function test_HeuristicCallback()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged = MOI.Bridges.Variable.Vectorize{Float64}(mock)
    attr = MOI.HeuristicCallback()
    f(::Any) = nothing
    MOI.set(bridged, attr, f)
    @test MOI.get(bridged, attr) === f
    return
end

function test_CallbackVariablePrimal()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged = MOI.Bridges.Variable.Vectorize{Float64}(mock)
    x, _ = MOI.add_constrained_variable(bridged, MOI.GreaterThan(1.0))
    y = MOI.get(mock, MOI.ListOfVariableIndices())[1]
    z = MOI.add_variable(bridged)
    attr = MOI.CallbackVariablePrimal(nothing)
    @test_throws(
        ErrorException("No mock callback primal is set for variable `$z`."),
        MOI.get(bridged, attr, z),
    )
    MOI.set(mock, attr, y, 1.0)
    MOI.set(mock, attr, z, 2.0)
    @test MOI.get(bridged, attr, z) == 2.0
    err = ArgumentError(
        "Variable bridge of type `$(typeof(MOI.Bridges.bridge(bridged, x)))` " *
        "does not support accessing the attribute `$attr`.",
    )
    @test_throws err MOI.get(bridged, attr, x)
    return
end

function test_LazyConstraint()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged = MOI.Bridges.Variable.Vectorize{Float64}(mock)
    x, _ = MOI.add_constrained_variable(bridged, MOI.GreaterThan(1.0))
    y = MOI.get(mock, MOI.ListOfVariableIndices())[1]
    sub = MOI.LazyConstraint(nothing)
    @test MOI.supports(bridged, sub)
    MOI.submit(bridged, sub, 2.0x, MOI.GreaterThan(1.0))
    mock.submitted[sub][1][1] ≈ 2.0y + 2.0
    mock.submitted[sub][1][2] == MOI.GreaterThan(1.0)
    return
end

function test_HeuristicSolution()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged = MOI.Bridges.Variable.Vectorize{Float64}(mock)
    x, _ = MOI.add_constrained_variable(bridged, MOI.GreaterThan(1.0))
    y = MOI.get(mock, MOI.ListOfVariableIndices())[1]
    z = MOI.add_variable(bridged)
    sub = MOI.HeuristicSolution(nothing)
    @test MOI.supports(bridged, sub)
    MOI.submit(bridged, sub, [z], [1.0])
    mock.submitted[sub][1][1] == [z]
    mock.submitted[sub][1][2] == [1.0]
    err = ErrorException(
        "Cannot substitute `$x` as it is bridged into `$(1.0y + 1.0)`.",
    )
    @test_throws err MOI.submit(bridged, sub, [x], [1.0])
    return
end

function test_NLPBlock()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged = MOI.Bridges.Variable.Vectorize{Float64}(mock)
    nlp_data =
        MOI.NLPBlockData([MOI.NLPBoundsPair(1.0, 2.0)], DummyEvaluator(), false)
    MOI.set(bridged, MOI.NLPBlock(), nlp_data)
    for nlp_data in
        [MOI.get(bridged, MOI.NLPBlock()), MOI.get(mock, MOI.NLPBlock())]
        @test nlp_data.constraint_bounds == [MOI.NLPBoundsPair(1.0, 2.0)]
        @test nlp_data.evaluator isa DummyEvaluator
        @test nlp_data.has_objective == false
    end
    return
end

# Model not supporting Interval
MOI.Utilities.@model(
    NoIntervalModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

function unsupported_constraint_attribute()
    mock = MOI.Utilities.MockOptimizer(NoIntervalModel{Float64}())
    bridged_mock = MOI.Bridges.Constraint.LessToGreater{Float64}(
        MOI.Bridges.Constraint.SplitInterval{Float64}(mock),
    )
    bridge = MOI.Bridges.Constraint.SplitIntervalBridge{
        Float64,
        MOI.VariableIndex,
        MOI.Interval{Float64},
        MOI.GreaterThan{Float64},
        MOI.LessThan{Float64},
    }
    attr = MOI.Test.UnknownConstraintAttribute()
    err = ArgumentError(
        "Bridge of type `$(bridge)` does not support accessing " *
        "the attribute `$attr`. If you encountered this error " *
        "unexpectedly, it probably means your model has been " *
        "reformulated using the bridge, and you are attempting to query " *
        "an attribute that we haven't implemented yet for this bridge. " *
        "Please open an issue at https://github.com/jump-dev/MathOptInterface.jl/issues/new " *
        "and provide a reproducible example explaining what you were " *
        "trying to do.",
    )
    x = MOI.add_variable(bridged_mock)
    ci = MOI.add_constraint(bridged_mock, x, MOI.Interval(0.0, 1.0))
    @test !MOI.Bridges.is_bridged(bridged_mock, ci)
    @test MOI.Bridges.is_bridged(bridged_mock.model, ci)
    @test !MOI.supports(bridged_mock, attr, typeof(ci))
    @test_throws err MOI.get(bridged_mock, attr, ci)
    return
end

function test_issue_453()
    mock = MOI.Utilities.MockOptimizer(NoIntervalModel{Float64}())
    bridged_mock = MOI.Bridges.Constraint.LessToGreater{Float64}(
        MOI.Bridges.Constraint.SplitInterval{Float64}(mock),
    )
    MOI.Utilities.loadfromstring!(
        bridged_mock,
        """
    variables: x
    maxobjective: 3.0x
    c: 2.0x in Interval(1.0, 4.0)
    x in LessThan(1.5)
""",
    )
    x = MOI.get(bridged_mock, MOI.VariableIndex, "x")
    @test isa(x, MOI.VariableIndex)
    c1 = MOI.get(
        bridged_mock,
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.Interval{Float64},
        },
        "c",
    )
    @test isa(
        c1,
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.Interval{Float64},
        },
    )
    c2 = MOI.get(bridged_mock, MOI.ConstraintIndex, "c")
    @test c1 == c2
    @test MOI.is_valid(
        bridged_mock,
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(x.value),
    )
    return
end

function test_custom_test()
    model = MOI.Bridges.Constraint.SplitInterval{Int}(NoIntervalModel{Int}())
    @test !MOI.Bridges.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{Float64},
        MOI.Interval{Float64},
    )

    x, y = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    f1 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 7)
    c1 = MOI.add_constraint(model, f1, MOI.Interval(-1, 1))

    @test MOI.get(model, MOI.ListOfConstraintTypesPresent()) ==
          [(MOI.ScalarAffineFunction{Int}, MOI.Interval{Int})]
    _test_num_constraints(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.GreaterThan{Int},
        0,
    )
    _test_num_constraints(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.Interval{Int},
        1,
    )
    @test (@inferred MOI.get(
        model,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Int},
            MOI.Interval{Int},
        }(),
    )) == [c1]

    f2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2, -1], [x, y]), 2)
    c2 = MOI.add_constraint(model, f1, MOI.GreaterThan(-2))

    @test MOI.get(model, MOI.ListOfConstraintTypesPresent()) == [
        (MOI.ScalarAffineFunction{Int}, MOI.GreaterThan{Int}),
        (MOI.ScalarAffineFunction{Int}, MOI.Interval{Int}),
    ]
    _test_num_constraints(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.GreaterThan{Int},
        1,
    )
    _test_num_constraints(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.Interval{Int},
        1,
    )
    @test (@inferred MOI.get(
        model,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Int},
            MOI.Interval{Int},
        }(),
    )) == [c1]
    @test (@inferred MOI.get(
        model,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Int},
            MOI.GreaterThan{Int},
        }(),
    )) == [c2]

    n = 4
    z = MOI.add_variables(model, n)
    scon_indices = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Int}}[]
    for (i, v) in enumerate([x; y; z])
        f = v
        c = MOI.add_constraint(model, f, MOI.Interval(i, 2i))
        push!(scon_indices, c)

        @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
            (MOI.ScalarAffineFunction{Int}, MOI.GreaterThan{Int}),
            (MOI.ScalarAffineFunction{Int}, MOI.Interval{Int}),
            (MOI.VariableIndex, MOI.Interval{Int}),
        ])
        _test_num_constraints(
            model,
            MOI.ScalarAffineFunction{Int},
            MOI.GreaterThan{Int},
            1,
        )
        _test_num_constraints(
            model,
            MOI.ScalarAffineFunction{Int},
            MOI.Interval{Int},
            1,
        )
        _test_num_constraints(model, MOI.VariableIndex, MOI.Interval{Int}, i)
        @test (@inferred MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Int},
                MOI.Interval{Int},
            }(),
        )) == [c1]
        @test (@inferred MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Int},
                MOI.GreaterThan{Int},
            }(),
        )) == [c2]
        # The indices should be returned in order of creation
        @test (@inferred MOI.get(
            model,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Interval{Int}}(),
        )) == scon_indices
    end

    vcon_indices = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}[]
    for (i, v) in enumerate(z)
        f = MOI.VectorOfVariables([v])
        c = MOI.add_constraint(model, f, MOI.Nonnegatives(1))
        push!(vcon_indices, c)

        @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
            (MOI.ScalarAffineFunction{Int}, MOI.GreaterThan{Int}),
            (MOI.ScalarAffineFunction{Int}, MOI.Interval{Int}),
            (MOI.VariableIndex, MOI.Interval{Int}),
            (MOI.VectorOfVariables, MOI.Nonnegatives),
        ])
        _test_num_constraints(
            model,
            MOI.ScalarAffineFunction{Int},
            MOI.GreaterThan{Int},
            1,
        )
        _test_num_constraints(
            model,
            MOI.ScalarAffineFunction{Int},
            MOI.Interval{Int},
            1,
        )
        _test_num_constraints(
            model,
            MOI.VariableIndex,
            MOI.Interval{Int},
            n + 2,
        )
        _test_num_constraints(model, MOI.VectorOfVariables, MOI.Nonnegatives, i)
        @test (@inferred MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Int},
                MOI.Interval{Int},
            }(),
        )) == [c1]
        @test (@inferred MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Int},
                MOI.GreaterThan{Int},
            }(),
        )) == [c2]
        # The indices should be returned in order of creation
        @test (@inferred MOI.get(
            model,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Interval{Int}}(),
        )) == scon_indices
        @test (@inferred MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonnegatives,
            }(),
        )) == vcon_indices
    end

    @test MOI.is_valid(model, c2)
    MOI.delete(model, c2)

    @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
        (MOI.ScalarAffineFunction{Int}, MOI.Interval{Int}),
        (MOI.VariableIndex, MOI.Interval{Int}),
        (MOI.VectorOfVariables, MOI.Nonnegatives),
    ])
    _test_num_constraints(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.GreaterThan{Int},
        0,
    )
    _test_num_constraints(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.Interval{Int},
        1,
    )
    _test_num_constraints(model, MOI.VariableIndex, MOI.Interval{Int}, n + 2)
    _test_num_constraints(model, MOI.VectorOfVariables, MOI.Nonnegatives, n)
    @test (@inferred MOI.get(
        model,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Int},
            MOI.Interval{Int},
        }(),
    )) == [c1]
    # The indices should be returned in order of creation
    @test (@inferred MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Interval{Int}}(),
    )) == scon_indices
    @test (@inferred MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonnegatives}(),
    )) == vcon_indices
    return
end

function test_MOI_Test()
    mock = MOI.Utilities.MockOptimizer(NoIntervalModel{Float64}())
    bridged_mock = MOI.Bridges.Constraint.LessToGreater{Float64}(
        MOI.Bridges.Constraint.SplitInterval{Float64}(mock),
    )
    MOI.Test.runtests(
        bridged_mock,
        MOI.Test.Config(exclude = Any[MOI.optimize!]),
        include = ["test_linear_", "test_model_"],
        exclude = [
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
            "test_model_ListOfConstraintAttributesSet",
        ],
    )
    return
end

function test_show()
    mock = MOI.Utilities.MockOptimizer(NoIntervalModel{Float64}())
    bridged_mock = MOI.Bridges.Constraint.LessToGreater{Float64}(
        MOI.Bridges.Constraint.SplitInterval{Float64}(mock),
    )
    @test_broken sprint(show, bridged_mock) == MOI.Utilities.replace_acronym(
        """
$(MOI.Bridges.Constraint.SingleBridgeOptimizer{MOI.Bridges.Constraint.LessToGreaterBridge{Float64,F,G} where G<:MOI.AbstractScalarFunction where F<:MOI.AbstractScalarFunction,MOI.Bridges.Constraint.SingleBridgeOptimizer{MOI.Bridges.Constraint.SplitIntervalBridge{Float64,F,S,LS,US} where US<:MOI.AbstractSet where LS<:MOI.AbstractSet where S<:MOI.AbstractSet where F<:MOI.AbstractFunction,MOI.Utilities.MockOptimizer{NoIntervalModel{Float64},Float64}}})
with 1 constraint bridge
with inner model $(MOI.Bridges.Constraint.SingleBridgeOptimizer{MOI.Bridges.Constraint.SplitIntervalBridge{Float64,F,S,LS,US} where US<:MOI.AbstractSet where LS<:MOI.AbstractSet where S<:MOI.AbstractSet where F<:MOI.AbstractFunction,MOI.Utilities.MockOptimizer{NoIntervalModel{Float64},Float64}})
  with 0 constraint bridges
  with inner model $(MOI.Utilities.MockOptimizer{NoIntervalModel{Float64},Float64})""",
    )
    return
end

MOI.Utilities.@model(
    AffineOnlyModel,
    (),
    (MOI.Interval,),
    (MOI.PositiveSemidefiniteConeTriangle,),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    (MOI.VectorAffineFunction,),
    true
)

function MOI.supports_constraint(
    ::AffineOnlyModel{T},
    ::Type{MOI.VariableIndex},
    ::Type{MOI.LessThan{T}},
) where {T}
    return false
end

function MOI.supports_constraint(
    ::AffineOnlyModel{T},
    ::Type{MOI.VariableIndex},
    ::Type{MOI.Interval{T}},
) where {T}
    return false
end

function test_double_deletion_scalar()
    # The variable is bridged to `VariableIndex`-in-`Interval` and then `ScalarAffineFunction`-in-`Interval`.
    # Hence there is two bridged `VariableIndex` constraints on the same variables and we need to be
    # careful not to delete the second one twice, see https://github.com/jump-dev/MathOptInterface.jl/issues/1231
    model =
        MOI.instantiate(AffineOnlyModel{Float64}, with_bridge_type = Float64)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, x, MOI.LessThan(1.0))
    # Need to test the bridging to make sure it's not functionized first as otherwise,
    # this test would not cover the case we want to test
    b1 = MOI.Bridges.bridge(model, c)
    @test b1 isa MOI.Bridges.Constraint.LessToIntervalBridge
    b2 = MOI.Bridges.bridge(model, b1.constraint)
    @test b2 isa MOI.Bridges.Constraint.ScalarFunctionizeBridge
    @test !MOI.Bridges.is_bridged(model, b2.constraint)
    MOI.delete(model, x)
    @test !MOI.is_valid(model, x)
    return
end

function test_double_deletion_vector()
    model =
        MOI.instantiate(AffineOnlyModel{Float64}, with_bridge_type = Float64)
    x = MOI.add_variables(model, 4)
    c = MOI.add_constraint(
        model,
        MOI.VectorOfVariables(x),
        MOI.PositiveSemidefiniteConeSquare(2),
    )
    b1 = MOI.Bridges.bridge(model, c)
    @test b1 isa MOI.Bridges.Constraint.SquareBridge
    b2 = MOI.Bridges.bridge(model, b1.triangle)
    @test b2 isa MOI.Bridges.Constraint.VectorFunctionizeBridge
    @test !MOI.Bridges.is_bridged(model, b2.constraint)
    MOI.delete(model, x)
    @test all(vi -> !MOI.is_valid(model, vi), x)
    return
end

function test_nesting_SingleBridgeOptimizer()
    T = Float64
    model = MOI.Utilities.Model{T}()
    b0 = MOI.Bridges.Variable.Free{T}(model)
    b1 = MOI.Bridges.Constraint.ScalarFunctionize{T}(b0)
    b2 = MOI.Bridges.Constraint.VectorFunctionize{T}(b1)
    b = MOI.Bridges.Objective.Functionize{T}(b2)
    x = MOI.add_variable(b)
    @test MOI.is_valid(b, x)
    @test !MOI.Bridges.is_bridged(b, x)
    @test MOI.is_valid(b2, x)
    @test !MOI.Bridges.is_bridged(b2, x)
    @test MOI.is_valid(b1, x)
    @test !MOI.Bridges.is_bridged(b1, x)
    @test MOI.is_valid(b0, x)
    @test MOI.Bridges.is_bridged(b0, x)
    @test !MOI.is_valid(model, x)
    clt = MOI.add_constraint(b, x, MOI.LessThan(one(T)))
    @test MOI.is_valid(b, clt)
    @test !MOI.Bridges.is_bridged(b, clt)
    @test MOI.is_valid(b2, clt)
    @test !MOI.Bridges.is_bridged(b2, clt)
    @test MOI.is_valid(b1, clt)
    @test MOI.Bridges.is_bridged(b1, clt)
    @test !MOI.is_valid(b0, clt)
    @test !MOI.is_valid(model, clt)
    cnn =
        MOI.add_constraint(b, MOI.Utilities.vectorize([x]), MOI.Nonnegatives(1))
    @test MOI.is_valid(b, cnn)
    @test !MOI.Bridges.is_bridged(b, cnn)
    @test MOI.is_valid(b2, cnn)
    @test MOI.Bridges.is_bridged(b2, cnn)
    @test !MOI.is_valid(b1, cnn)
    @test !MOI.is_valid(b0, cnn)
    @test !MOI.is_valid(model, cnn)
    obj = x
    attr = MOI.ObjectiveFunction{typeof(obj)}()
    MOI.set(b, attr, obj)
    @test MOI.Bridges.is_bridged(b, attr)
    obj = one(T) * x
    attr = MOI.ObjectiveFunction{typeof(obj)}()
    MOI.set(b, attr, obj)
    @test !MOI.Bridges.is_bridged(b, attr)
    @test !MOI.Bridges.is_bridged(b2, attr)
    @test !MOI.Bridges.is_bridged(b1, attr)
    @test !MOI.Bridges.is_bridged(b0, attr)
    return
end

struct UnsupportedSet <: MOI.AbstractScalarSet end

function test_supports_bridging_constrained_variable()
    model = MOI.Bridges.Constraint.SplitInterval{Int}(NoIntervalModel{Int}())
    @test !MOI.Bridges.supports_bridging_constrained_variable(
        model,
        UnsupportedSet,
    )
    return
end

function test_supports_add_constrained_variable()
    model = MOI.Bridges.Constraint.SplitInterval{Int}(NoIntervalModel{Int}())
    @test !MOI.supports_add_constrained_variable(model, UnsupportedSet)
    return
end

function test_supports_bridging_constraint()
    model = MOI.Bridges.Constraint.SplitInterval{Int}(NoIntervalModel{Int}())
    @test !MOI.Bridges.supports_bridging_constraint(
        model,
        MOI.VariableIndex,
        UnsupportedSet,
    )
    return
end

function test_supports_bridging_objective_function()
    model = MOI.Bridges.Constraint.SplitInterval{Int}(NoIntervalModel{Int}())
    @test !MOI.Bridges.supports_bridging_objective_function(
        model,
        MOI.VariableIndex,
    )
    return
end

function test_get_ObjectiveFunctionType()
    model = MOI.Bridges.Constraint.SplitInterval{Int}(NoIntervalModel{Int}())
    x = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    @test MOI.get(model, MOI.ObjectiveFunctionType()) == MOI.VariableIndex
    return
end

include("identity_bridge.jl")

function test_recursive_model_variable(::Type{T} = Int) where {T}
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}())
    BT = IdentityBridges.VariableBridge{T}
    b = MOI.Bridges.Variable.SingleBridgeOptimizer{BT}(model)
    x, cx = MOI.add_constrained_variable(b, MOI.EqualTo(one(T)))
    @test MOI.Bridges.is_bridged(b, x)
    @test MOI.Bridges.is_bridged(b, cx)
    @test MOI.get(b, MOI.ConstraintFunction(), cx) == x
    @test MOI.get(b, MOI.ConstraintSet(), cx) == MOI.EqualTo(one(T))
    MOI.set(b, MOI.ConstraintSet(), cx, MOI.EqualTo(zero(T)))
    @test MOI.get(b, MOI.ConstraintSet(), cx) == MOI.EqualTo(zero(T))
    @test MOI.is_valid(b, x)
    MOI.delete(b, x)
    @test !MOI.is_valid(b, x)
end

function test_recursive_model_constraint(::Type{T} = Int) where {T}
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}())
    BT = IdentityBridges.ConstraintBridge{T}
    b = MOI.Bridges.Constraint.SingleBridgeOptimizer{BT}(model)
    x = MOI.add_variable(b)
    func = one(T) * x
    set = MOI.EqualTo(zero(T))
    c = MOI.add_constraint(b, func, set)
    @test MOI.Bridges.is_bridged(b, c)
    @test MOI.get(b, MOI.ConstraintFunction(), c) ≈ func
    new_func = T(2) * x
    MOI.set(b, MOI.ConstraintFunction(), c, new_func)
    @test MOI.get(b, MOI.ConstraintFunction(), c) == new_func
    MOI.modify(b, c, MOI.ScalarCoefficientChange(x, T(3)))
    @test MOI.get(b, MOI.ConstraintFunction(), c) ≈ T(3) * x
    MOI.modify(b, c, MOI.ScalarConstantChange(T(-1)))
    @test MOI.get(b, MOI.ConstraintFunction(), c) ≈ T(3) * x + T(-1)
    @test MOI.get(b, MOI.ConstraintSet(), c) == set
    new_set = MOI.EqualTo(one(T))
    MOI.set(b, MOI.ConstraintSet(), c, new_set)
    @test MOI.get(b, MOI.ConstraintSet(), c) == new_set
    MOI.set(b, MOI.ConstraintDualStart(), c, one(T))
    @test MOI.get(b, MOI.ConstraintDualStart(), c) == one(T)
    @test MOI.is_valid(b, c)
    MOI.delete(b, c)
    @test !MOI.is_valid(b, c)
end

function test_recursive_model_objective(::Type{T} = Int) where {T}
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}())
    BT = IdentityBridges.ObjectiveBridge{T}
    b = MOI.Bridges.Objective.SingleBridgeOptimizer{BT}(model)
    x = MOI.add_variable(b)
    func = one(T) * x
    MOI.set(b, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test !MOI.Bridges.is_objective_bridged(b)
    attr = MOI.ObjectiveFunction{typeof(func)}()
    MOI.set(b, attr, func)
    @test MOI.Bridges.is_objective_bridged(b)
    @test MOI.get(b, attr) ≈ func
    attr = MOI.ObjectiveFunction{typeof(x)}()
    MOI.set(b, attr, x)
    @test !MOI.Bridges.is_objective_bridged(b)
    @test MOI.get(b, attr) ≈ x
end

function test_invalid_modifications()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.Model{Float64}(),
        Float64,
    )
    config = MOI.Test.Config()
    MOI.Test.test_modification_set_function_single_variable(model, config)
    MOI.empty!(model)
    MOI.Test.test_modification_incorrect(model, config)
    MOI.empty!(model)
    MOI.Test.test_modification_incorrect_VariableIndex(model, config)
    return
end

function test_ListOfConstraintAttributesSet()
    model = MOI.Bridges.Constraint.SplitInterval{Float64}(
        MOI.Utilities.Model{Float64}(),
    )
    x = MOI.add_variable(model)
    f = 1.0 * x
    s = MOI.Interval(0.0, 1.0)
    c = MOI.add_constraint(model, f, s)
    MOI.set(model, MOI.ConstraintName(), c, "c")
    @test MOI.get(
        model,
        MOI.ListOfConstraintAttributesSet{typeof(f),typeof(s)}(),
    ) == [MOI.ConstraintName()]
    return
end

end  # module

TestBridgeOptimizer.runtests()


text:
"""
    binSearch(array, x)

Find the smallest index `i` in sorted vector `array` that satisfies `array[i] <= x < array[i+1]`.

# Examples

```julia-repl
julia> binSearch([1,2.5,3.4,4,5], 3.2)
2
```
"""
function binSearch(array::Vector{T}, x::U)::Int where {T,U <: Real}
    pos_start = 1; pos_end = length(array);

    while pos_start < pos_end
        if pos_start == pos_end - 1
            if x < array[pos_end]
                return pos_start
            end
        end

        pos_middle = (pos_start + pos_end)÷2
        if x < array[pos_middle]
            pos_end = pos_middle
        else
            pos_start = pos_middle
        end
    end

    pos_start
end

ext: jl


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    AbstractNeighborhood

A neighborhood on a spatial domain of type `D`.
"""
abstract type AbstractNeighborhood{D<:AbstractDomain} end

# Neighborhoods are functor objects that can be evaluated
# at a given location:
#
# julia> neighborhood(location)
#
# The operator () returns the neighbors (as integers) and
# is implemented differently depending on the domain type.

#------------------
# IMPLEMENTATIONS
#------------------
include("neighborhoods/cube_neighborhood.jl")
include("neighborhoods/ball_neighborhood.jl")


text:
# Libraries
using PGFPlotsX, LaTeXStrings;
using MessyTimeSeries: no_combinations;
using Statistics;

# Custom functions for AJK and BJK probabilities

# - AJK
compute_ajk_cases(nT::Int64, d::Int64, k::Int64) = no_combinations(d, k)*no_combinations(nT-d, d-k);
compute_ajk_probability(nT::Int64, d::Int64, k::Int64) = compute_ajk_cases(nT, d, k)/no_combinations(nT, d);

# - BJK
compute_bjk_ith_cases(i::Int64, no_subsamples::Int64, c::Int64, k::Int64) = (no_subsamples-i >= c-k) + (i-1 >= c-k);
compute_bjk_probability(no_subsamples::Int64, c::Int64, k::Int64) = mapreduce(i -> compute_bjk_ith_cases(i, no_subsamples, c, k), +, 1:no_subsamples)/(no_subsamples^2);

#=
Custom functions to compute the expected value of the Jaccard indices
=#

# - Main methods
expected_jaccard_ajk(nT::Int64, d::Int64) = mapreduce(k -> k/(2*d-k)*compute_ajk_probability(nT, d, k), +, 1:d-1);
expected_jaccard_bjk(T::Int64, c::Int64)  = mapreduce(k -> k/(2*c-k)*compute_bjk_probability(T-c+1, c, k), +, 1:c-1);

# - Promotion rules
expected_jaccard_ajk(nT::Real, d::Real) = expected_jaccard_ajk(Int64(nT), Int64(floor(d)));
expected_jaccard_bjk(T::Real, c::Real) = expected_jaccard_bjk(Int64(T), Int64(floor(c)));

#=
Surface chart
=#

# Plot: grid
push!(PGFPlotsX.CUSTOM_PREAMBLE,
        raw"\usepgfplotslibrary{colorbrewer}",
        raw"\usepgfplotslibrary{colormaps}",
        raw"\usepgfplotslibrary{patchplots}",
        raw"\pgfplotsset
            {
                tick label style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
                label style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
                legend style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
                title style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
            }"
)

# GroupPlot axs
gpr = Array{Any,1}(undef, 5);

# Loop over c
for (i, c) in enumerate([4, 8, 12, 16, 20])

    # n,T grids
    n_grid = range(5, stop=50, length=10);
    T_grid = range(50, stop=500, length=10);

    # Compute z_grid
    title_string = L"c=%$(c)";
    z_grid = [expected_jaccard_bjk(T, c)-expected_jaccard_ajk(n*T, n*c) for n in n_grid, T in T_grid];

    # Plots
    p1 = @pgf Plot3({surf, point_meta="explicit"}, Coordinates(n_grid, T_grid, z_grid, meta=z_grid));

    gpr[i] = @pgf Axis(
    {
        "axis on top"="false",
        height = "10cm",
        width = "10cm",
        grid = "major",
        grid_style = "{darkgray, very thin}",
        tick_style = "{darkgray, very thin}",
        "colormap"="{turbo}{rgb=(0.18995,0.07176,0.23217), rgb=(0.19483,0.08339,0.26149), rgb=(0.19956,0.09498,0.29024), rgb=(0.20415,0.10652,0.31844), rgb=(0.20860,0.11802,0.34607), rgb=(0.21291,0.12947,0.37314), rgb=(0.21708,0.14087,0.39964), rgb=(0.22111,0.15223,0.42558), rgb=(0.22500,0.16354,0.45096), rgb=(0.22875,0.17481,0.47578), rgb=(0.23236,0.18603,0.50004), rgb=(0.23582,0.19720,0.52373), rgb=(0.23915,0.20833,0.54686), rgb=(0.24234,0.21941,0.56942), rgb=(0.24539,0.23044,0.59142), rgb=(0.24830,0.24143,0.61286), rgb=(0.25107,0.25237,0.63374), rgb=(0.25369,0.26327,0.65406), rgb=(0.25618,0.27412,0.67381), rgb=(0.25853,0.28492,0.69300), rgb=(0.26074,0.29568,0.71162), rgb=(0.26280,0.30639,0.72968), rgb=(0.26473,0.31706,0.74718), rgb=(0.26652,0.32768,0.76412), rgb=(0.26816,0.33825,0.78050), rgb=(0.26967,0.34878,0.79631), rgb=(0.27103,0.35926,0.81156), rgb=(0.27226,0.36970,0.82624), rgb=(0.27334,0.38008,0.84037), rgb=(0.27429,0.39043,0.85393), rgb=(0.27509,0.40072,0.86692), rgb=(0.27576,0.41097,0.87936), rgb=(0.27628,0.42118,0.89123), rgb=(0.27667,0.43134,0.90254), rgb=(0.27691,0.44145,0.91328), rgb=(0.27701,0.45152,0.92347), rgb=(0.27698,0.46153,0.93309), rgb=(0.27680,0.47151,0.94214), rgb=(0.27648,0.48144,0.95064), rgb=(0.27603,0.49132,0.95857), rgb=(0.27543,0.50115,0.96594), rgb=(0.27469,0.51094,0.97275), rgb=(0.27381,0.52069,0.97899), rgb=(0.27273,0.53040,0.98461), rgb=(0.27106,0.54015,0.98930), rgb=(0.26878,0.54995,0.99303), rgb=(0.26592,0.55979,0.99583), rgb=(0.26252,0.56967,0.99773), rgb=(0.25862,0.57958,0.99876), rgb=(0.25425,0.58950,0.99896), rgb=(0.24946,0.59943,0.99835), rgb=(0.24427,0.60937,0.99697), rgb=(0.23874,0.61931,0.99485), rgb=(0.23288,0.62923,0.99202), rgb=(0.22676,0.63913,0.98851), rgb=(0.22039,0.64901,0.98436), rgb=(0.21382,0.65886,0.97959), rgb=(0.20708,0.66866,0.97423), rgb=(0.20021,0.67842,0.96833), rgb=(0.19326,0.68812,0.96190), rgb=(0.18625,0.69775,0.95498), rgb=(0.17923,0.70732,0.94761), rgb=(0.17223,0.71680,0.93981), rgb=(0.16529,0.72620,0.93161), rgb=(0.15844,0.73551,0.92305), rgb=(0.15173,0.74472,0.91416), rgb=(0.14519,0.75381,0.90496), rgb=(0.13886,0.76279,0.89550), rgb=(0.13278,0.77165,0.88580), rgb=(0.12698,0.78037,0.87590), rgb=(0.12151,0.78896,0.86581), rgb=(0.11639,0.79740,0.85559), rgb=(0.11167,0.80569,0.84525), rgb=(0.10738,0.81381,0.83484), rgb=(0.10357,0.82177,0.82437), rgb=(0.10026,0.82955,0.81389), rgb=(0.09750,0.83714,0.80342), rgb=(0.09532,0.84455,0.79299), rgb=(0.09377,0.85175,0.78264), rgb=(0.09287,0.85875,0.77240), rgb=(0.09267,0.86554,0.76230), rgb=(0.09320,0.87211,0.75237), rgb=(0.09451,0.87844,0.74265), rgb=(0.09662,0.88454,0.73316), rgb=(0.09958,0.89040,0.72393), rgb=(0.10342,0.89600,0.71500), rgb=(0.10815,0.90142,0.70599), rgb=(0.11374,0.90673,0.69651), rgb=(0.12014,0.91193,0.68660), rgb=(0.12733,0.91701,0.67627), rgb=(0.13526,0.92197,0.66556), rgb=(0.14391,0.92680,0.65448), rgb=(0.15323,0.93151,0.64308), rgb=(0.16319,0.93609,0.63137), rgb=(0.17377,0.94053,0.61938), rgb=(0.18491,0.94484,0.60713), rgb=(0.19659,0.94901,0.59466), rgb=(0.20877,0.95304,0.58199), rgb=(0.22142,0.95692,0.56914), rgb=(0.23449,0.96065,0.55614), rgb=(0.24797,0.96423,0.54303), rgb=(0.26180,0.96765,0.52981), rgb=(0.27597,0.97092,0.51653), rgb=(0.29042,0.97403,0.50321), rgb=(0.30513,0.97697,0.48987), rgb=(0.32006,0.97974,0.47654), rgb=(0.33517,0.98234,0.46325), rgb=(0.35043,0.98477,0.45002), rgb=(0.36581,0.98702,0.43688), rgb=(0.38127,0.98909,0.42386), rgb=(0.39678,0.99098,0.41098), rgb=(0.41229,0.99268,0.39826), rgb=(0.42778,0.99419,0.38575), rgb=(0.44321,0.99551,0.37345), rgb=(0.45854,0.99663,0.36140), rgb=(0.47375,0.99755,0.34963), rgb=(0.48879,0.99828,0.33816), rgb=(0.50362,0.99879,0.32701), rgb=(0.51822,0.99910,0.31622), rgb=(0.53255,0.99919,0.30581), rgb=(0.54658,0.99907,0.29581), rgb=(0.56026,0.99873,0.28623), rgb=(0.57357,0.99817,0.27712), rgb=(0.58646,0.99739,0.26849), rgb=(0.59891,0.99638,0.26038), rgb=(0.61088,0.99514,0.25280), rgb=(0.62233,0.99366,0.24579), rgb=(0.63323,0.99195,0.23937), rgb=(0.64362,0.98999,0.23356), rgb=(0.65394,0.98775,0.22835), rgb=(0.66428,0.98524,0.22370), rgb=(0.67462,0.98246,0.21960), rgb=(0.68494,0.97941,0.21602), rgb=(0.69525,0.97610,0.21294), rgb=(0.70553,0.97255,0.21032), rgb=(0.71577,0.96875,0.20815), rgb=(0.72596,0.96470,0.20640), rgb=(0.73610,0.96043,0.20504), rgb=(0.74617,0.95593,0.20406), rgb=(0.75617,0.95121,0.20343), rgb=(0.76608,0.94627,0.20311), rgb=(0.77591,0.94113,0.20310), rgb=(0.78563,0.93579,0.20336), rgb=(0.79524,0.93025,0.20386), rgb=(0.80473,0.92452,0.20459), rgb=(0.81410,0.91861,0.20552), rgb=(0.82333,0.91253,0.20663), rgb=(0.83241,0.90627,0.20788), rgb=(0.84133,0.89986,0.20926), rgb=(0.85010,0.89328,0.21074), rgb=(0.85868,0.88655,0.21230), rgb=(0.86709,0.87968,0.21391), rgb=(0.87530,0.87267,0.21555), rgb=(0.88331,0.86553,0.21719), rgb=(0.89112,0.85826,0.21880), rgb=(0.89870,0.85087,0.22038), rgb=(0.90605,0.84337,0.22188), rgb=(0.91317,0.83576,0.22328), rgb=(0.92004,0.82806,0.22456), rgb=(0.92666,0.82025,0.22570), rgb=(0.93301,0.81236,0.22667), rgb=(0.93909,0.80439,0.22744), rgb=(0.94489,0.79634,0.22800), rgb=(0.95039,0.78823,0.22831), rgb=(0.95560,0.78005,0.22836), rgb=(0.96049,0.77181,0.22811), rgb=(0.96507,0.76352,0.22754), rgb=(0.96931,0.75519,0.22663), rgb=(0.97323,0.74682,0.22536), rgb=(0.97679,0.73842,0.22369), rgb=(0.98000,0.73000,0.22161), rgb=(0.98289,0.72140,0.21918), rgb=(0.98549,0.71250,0.21650), rgb=(0.98781,0.70330,0.21358), rgb=(0.98986,0.69382,0.21043), rgb=(0.99163,0.68408,0.20706), rgb=(0.99314,0.67408,0.20348), rgb=(0.99438,0.66386,0.19971), rgb=(0.99535,0.65341,0.19577), rgb=(0.99607,0.64277,0.19165), rgb=(0.99654,0.63193,0.18738), rgb=(0.99675,0.62093,0.18297), rgb=(0.99672,0.60977,0.17842), rgb=(0.99644,0.59846,0.17376), rgb=(0.99593,0.58703,0.16899), rgb=(0.99517,0.57549,0.16412), rgb=(0.99419,0.56386,0.15918), rgb=(0.99297,0.55214,0.15417), rgb=(0.99153,0.54036,0.14910), rgb=(0.98987,0.52854,0.14398), rgb=(0.98799,0.51667,0.13883), rgb=(0.98590,0.50479,0.13367), rgb=(0.98360,0.49291,0.12849), rgb=(0.98108,0.48104,0.12332), rgb=(0.97837,0.46920,0.11817), rgb=(0.97545,0.45740,0.11305), rgb=(0.97234,0.44565,0.10797), rgb=(0.96904,0.43399,0.10294), rgb=(0.96555,0.42241,0.09798), rgb=(0.96187,0.41093,0.09310), rgb=(0.95801,0.39958,0.08831), rgb=(0.95398,0.38836,0.08362), rgb=(0.94977,0.37729,0.07905), rgb=(0.94538,0.36638,0.07461), rgb=(0.94084,0.35566,0.07031), rgb=(0.93612,0.34513,0.06616), rgb=(0.93125,0.33482,0.06218), rgb=(0.92623,0.32473,0.05837), rgb=(0.92105,0.31489,0.05475), rgb=(0.91572,0.30530,0.05134), rgb=(0.91024,0.29599,0.04814), rgb=(0.90463,0.28696,0.04516), rgb=(0.89888,0.27824,0.04243), rgb=(0.89298,0.26981,0.03993), rgb=(0.88691,0.26152,0.03753), rgb=(0.88066,0.25334,0.03521), rgb=(0.87422,0.24526,0.03297), rgb=(0.86760,0.23730,0.03082), rgb=(0.86079,0.22945,0.02875), rgb=(0.85380,0.22170,0.02677), rgb=(0.84662,0.21407,0.02487), rgb=(0.83926,0.20654,0.02305), rgb=(0.83172,0.19912,0.02131), rgb=(0.82399,0.19182,0.01966), rgb=(0.81608,0.18462,0.01809), rgb=(0.80799,0.17753,0.01660), rgb=(0.79971,0.17055,0.01520), rgb=(0.79125,0.16368,0.01387), rgb=(0.78260,0.15693,0.01264), rgb=(0.77377,0.15028,0.01148), rgb=(0.76476,0.14374,0.01041), rgb=(0.75556,0.13731,0.00942), rgb=(0.74617,0.13098,0.00851), rgb=(0.73661,0.12477,0.00769), rgb=(0.72686,0.11867,0.00695), rgb=(0.71692,0.11268,0.00629), rgb=(0.70680,0.10680,0.00571), rgb=(0.69650,0.10102,0.00522), rgb=(0.68602,0.09536,0.00481), rgb=(0.67535,0.08980,0.00449), rgb=(0.66449,0.08436,0.00424), rgb=(0.65345,0.07902,0.00408), rgb=(0.64223,0.07380,0.00401), rgb=(0.63082,0.06868,0.00401), rgb=(0.61923,0.06367,0.00410), rgb=(0.60746,0.05878,0.00427), rgb=(0.59550,0.05399,0.00453), rgb=(0.58336,0.04931,0.00486), rgb=(0.57103,0.04474,0.00529), rgb=(0.55852,0.04028,0.00579), rgb=(0.54583,0.03593,0.00638), rgb=(0.53295,0.03169,0.00705), rgb=(0.51989,0.02756,0.00780), rgb=(0.50664,0.02354,0.00863), rgb=(0.49321,0.01963,0.00955), rgb=(0.47960,0.01583,0.01055)}",
        "colorbar"=ifelse(i==5, "true", "false"),
        "colorbar style"=raw"
        {
            ytick={0,0.02,...,0.12},
            yticklabel style=
            {
                text width=2.5em,
                align=right,
                /pgf/number format/.cd,
                    fixed,
                    fixed zerofill
            }
        }",
        "colorbar/width" = "0.8cm",
        title = title_string,
        view = (45, 15),
        point_meta_min = "0.0",
        point_meta_max = "0.12",
        zmin = "0",
        zmax = "0.15",
        "scaled z ticks"="false",
        "zticklabel"=raw"\pgfkeys{/pgf/number format/.cd,fixed,precision=2,zerofill}\pgfmathprintnumber{\tick}",
        xlabel = L"n",
        ylabel = L"T",
        zlabel = ifelse(i==1, "Expected Jaccard Index (difference)", ""),
    },
        p1,
    );
end

fig = @pgf TikzPicture(GroupPlot(
    { group_style = {group_size="5 by 1", vertical_sep="50pt", horizontal_sep="50pt"},
      height="200pt",
      width="225pt"
    },
    gpr...,)
)

pgfsave("./img/appendix/simulations_surface.pdf", fig);


text:
import ..TestSuite

@testset "Unit normalization" begin
    TestSuite.runtests(LegendreUnitNorm())
end
@testset "4π normalization" begin
    TestSuite.runtests(LegendreFourPiNorm())
end
@testset "Orthonormal normalization" begin
    TestSuite.runtests(LegendreOrthoNorm())
end
@testset "Spherical normalization" begin
    TestSuite.runtests(LegendreSphereNorm())
end


text:
using Parameters
# using PyCall
using PyPlot
# using VideoIO
include("bilog.jl")
include("colorize_field.jl")

# anim = pyimport("matplotlib.animation")  # After using PyPlot, so matplotlib installed via conda.
# const wave_propagation_filename = "wave_propagation.mp4"
const beam_energy_map_filename = "beam_energy_map.png"
const transmit_time_map_filename = "transmit_time_map.png"

function saveall(images, beam_energy_map, transmit_time_map, sim_params, output_path="images")
    @unpack fov, dbrange, orientation = sim_params

    extent = [0, fov[2], -1/2 * fov[1], 1/2 * fov[1]]
    xlabel = "Depth [m]"
    ylabel = "Azimuth [m]"

    vertical_enabled = orientation == :vertical
    if vertical_enabled
        extent = [-1/2 * fov[1], 1/2 * fov[1], fov[2], 0]
        xlabel, ylabel = ylabel, xlabel

        images = mapslices(rotr90, images; dims=[1, 2])
        beam_energy_map = rotr90(beam_energy_map)
        transmit_time_map = rotr90(transmit_time_map)
    end

    # wave_field = bilog(images, dbrange)
    # vmin, vmax = extrema(wave_field)
    #
    # fig = figure()
    # ims = []
    # for i_time in 1:size(wave_field, 3)
    #     im = PyPlot.imshow(wave_field[:, :, i_time], extent=extent, cmap="viridis", vmin=vmin, vmax=vmax)
    #     if 1 == i_time
    #         PyPlot.colorbar()
    #     end
    #     PyPlot.title("Wave amplitude [dB]")
    #     PyPlot.xlabel(xlabel)
    #     PyPlot.ylabel(ylabel)
    #     push!(ims, PyCall.PyObject[im])
    # end
    # #= close() =#
    #
    # # If matplotlib complains, ensure that
    # # a) ffmpeg is installed with libx264 support, and
    # # b) matplotlib is built with ffmpeg support enabled
    # ani = anim.ArtistAnimation(fig, ims, interval=30, blit=true, repeat=false)
    # ffmpeg_writer = anim.FFMpegWriter()
    # ani.save(joinpath(output_path, wave_propagation_filename), ffmpeg_writer, dpi=72)
    # # ani.save("ani.gif");

    PyPlot.figure()
    PyPlot.imshow(bilog(beam_energy_map, dbrange), extent=extent)
    PyPlot.colorbar()
    PyPlot.title("Beam energy map [dB]")
    PyPlot.xlabel(xlabel)
    PyPlot.ylabel(ylabel)
    PyPlot.savefig(joinpath(output_path, beam_energy_map_filename))
    PyPlot.close()

    PyPlot.figure()
    PyPlot.imshow(transmit_time_map .* 1e6, extent=extent)
    PyPlot.colorbar()
    PyPlot.title("Transmit time map [µs]")
    PyPlot.xlabel(xlabel)
    PyPlot.ylabel(ylabel)
    PyPlot.savefig(joinpath(output_path, transmit_time_map_filename))
    PyPlot.close()

    return  # nothing
end


text:
# ---
# title: 1665. Minimum Initial Energy to Finish Tasks
# id: problem1665
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-initial-energy-to-finish-tasks/description/>
# hidden: true
# ---
#
# You are given an array `tasks` where `tasks[i] = [actuali, minimumi]`:
#
#   * `actuali` is the actual amount of energy you **spend to finish** the `ith` task.
#   * `minimumi` is the minimum amount of energy you **require to begin** the `ith` task.
#
# For example, if the task is `[10, 12]` and your current energy is `11`, you
# cannot start this task. However, if your current energy is `13`, you can
# complete this task, and your energy will be `3` after finishing it.
#
# You can finish the tasks in **any order** you like.
#
# Return _the **minimum** initial amount of energy you will need_ _to finish all
# the tasks_.
#
#
#
# **Example 1:**
#
#
#
#     Input: tasks = [[1,2],[2,4],[4,8]]
#     Output: 8
#     Explanation:
#     Starting with 8 energy, we finish the tasks in the following order:
#         - 3rd task. Now energy = 8 - 4 = 4.
#         - 2nd task. Now energy = 4 - 2 = 2.
#         - 1st task. Now energy = 2 - 1 = 1.
#     Notice that even though we have leftover energy, starting with 7 energy does not work because we cannot do the 3rd task.
#
# **Example 2:**
#
#
#
#     Input: tasks = [[1,3],[2,4],[10,11],[10,12],[8,9]]
#     Output: 32
#     Explanation:
#     Starting with 32 energy, we finish the tasks in the following order:
#         - 1st task. Now energy = 32 - 1 = 31.
#         - 2nd task. Now energy = 31 - 2 = 29.
#         - 3rd task. Now energy = 29 - 10 = 19.
#         - 4th task. Now energy = 19 - 10 = 9.
#         - 5th task. Now energy = 9 - 8 = 1.
#
# **Example 3:**
#
#
#
#     Input: tasks = [[1,7],[2,8],[3,9],[4,10],[5,11],[6,12]]
#     Output: 27
#     Explanation:
#     Starting with 27 energy, we finish the tasks in the following order:
#         - 5th task. Now energy = 27 - 5 = 22.
#         - 2nd task. Now energy = 22 - 2 = 20.
#         - 3rd task. Now energy = 20 - 3 = 17.
#         - 1st task. Now energy = 17 - 1 = 16.
#         - 4th task. Now energy = 16 - 4 = 12.
#         - 6th task. Now energy = 12 - 6 = 6.
#
#
#
#
# **Constraints:**
#
#   * `1 <= tasks.length <= 105`
#   * `1 <= actuali <= minimumi <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
mutable struct TemporalGraph
    #graph::SimpleWeightedDiGraph
    graph::Dict
    nv::Int64
    steps::Int64
    # Visual coordinates
    locs_x::Vector{Float64}
    locs_y::Vector{Float64}
    has_async_nodes::Bool

    TemporalGraph() = new(Dict{String,SimpleWeightedDiGraph}(), 0, 1, Vector{Float64}(), Vector{Float64}(), false)
end

function TemporalGraph(network::QNetwork, steps::Int64; memory_prob::Float64=1.0,
    memory_costs::Dict=Dict())::TemporalGraph
    temp_graph = TemporalGraph()
    temp_graph.nv = length(network.nodes)
    temp_graph.steps = steps
    temp_graph.has_async_nodes = false

    # Make a copy of the network so we don't change QNetwork structure.
    netcopy = deepcopy(network)

    if memory_prob != 1.0
        @assert 0 <= memory_prob <= 1
        # Iterate through nodes and randomly reassign memories according to memory_prob
        # TODO: BAD PRACTICE. Don't manipulate QNetwork data.
        for node in netcopy.nodes
            if rand(Float64) <= memory_prob
                node.has_memory = true
            else
                node.has_memory = false
            end
        end
    end

    # Keep tabs on which nodes have memories and what the the associated cost are
    node_memories = Dict()
    node_memory_costs = Dict()
    for node in netcopy.nodes
        node_memories[node.id] = node.has_memory
        node_memory_costs[node.id] = node.memory_costs
    end

    for cost_key in keys(zero_costvector())
        temp_graph.graph[cost_key] = SimpleWeightedDiGraph()

        for t in 1:steps
            # Vertices
            add_vertices!(temp_graph.graph[cost_key], temp_graph.nv)

            # Channels
            for channel in netcopy.channels
                if channel.active == true
                    src = findfirst(x -> x == channel.src, netcopy.nodes) + (t-1)*temp_graph.nv
                    dest = findfirst(x -> x == channel.dest, netcopy.nodes) + (t-1)*temp_graph.nv
                    weight = channel.costs[cost_key]
                    add_edge!(temp_graph.graph[cost_key], src, dest, weight)
                    add_edge!(temp_graph.graph[cost_key], dest, src, weight)
                end
            end
        end

        # Memory channels
        for t in 1:(steps-1)
            for node in 1:temp_graph.nv
                # Add temporal link if memory exists
                if node_memories[node] == true
                    # WARNING: Memory costs aren't dependent on the timestep size.
                    # If the user specified a memory cost, use that.
                    if length(memory_costs) != 0
                        cost = memory_costs[cost_key]
                    # Otherwise pull it from the node attribute
                    else
                        cv = node_memory_costs[node]
                        cost = cv[cost_key]
                    end
                    # If the cost is zero, set it to epsilon so SimpleWeightedGraphs
                    # add it to the sparse adj. matrix.
                    if cost == 0
                        cost = eps(Float64)
                    end

                    src = node + (t-1)*temp_graph.nv
                    dest = src + temp_graph.nv
                    add_edge!(temp_graph.graph[cost_key], src, dest, cost)
                end
            end
        end

        # Coords

        # Offset factor for temporal nodes
        offsetX = 0.1
        offsetY = 0.2

        for t in 1:steps
            for node in netcopy.nodes
                push!(temp_graph.locs_x, node.location.x + (t-1) * offsetX)
                push!(temp_graph.locs_y, node.location.y + (t-1) * offsetY)
            end
        end
    end

    return temp_graph
end

"""
Add asynchronus nodes to a Temporal Network.

This function adds tempnet.nv nodes to the graph. If a given node is source,
outgoing asynchronus edges are added to the node that connects it to its temporal
counterparts. This is done with incremental timeweights ϵ so that earlier times
are prioritized.

Likewise if a node is a dest, incoming asynchronus edges are added with incremental
weights ϵ.
"""
function add_async_nodes!(tempnet::QuNet.TemporalGraph, endusers::Vector{Tuple{Int64, Int64}};
    ϵ=1)

    # No double dipping!
    if tempnet.has_async_nodes == true
        @warn("Redundent use of add_async_nodes -- T.has_async_nodes == true
        Returning an unmodified TemporalGraph")
        return tempnet
    end

    N = tempnet.nv
    steps = tempnet.steps
    # index offset for asynchronus nodes
    async_offset = N*steps

    # Keep track of asynchronus src and dst nodes
    src_nodes = []
    dst_nodes = []
    for userpair in endusers
        src = userpair[1]
        dst = userpair[2]
        if src > async_offset
            push!(src_nodes, src)
        end
        if dst > async_offset
            push!(dst_nodes, dst)
        end
    end

    for costkey in keys(zero_costvector())
        graph = tempnet.graph[costkey]
        # Add the asynchronus nodes
        add_vertices!(graph, N)

        # Connect async nodes to temporal counterparts:
        for async in (1 + async_offset):(N + async_offset)
            # If node is a src, make asynchronus channels outgoing
            if async in src_nodes
                for t in 1:steps
                    # Connect each temporal node to its async counterpart
                    # Use time_dependent weight to prioritise top layers when routing
                    #add_edge(src, dst)
                    add_edge!(graph, async, async%async_offset + N*(t-1), ϵ*t)
                end
            # If node is a dst, make asynchronus channels incoming
            elseif async in dst_nodes
                for t in 1:steps
                    add_edge!(graph, async%async_offset + N*(t-1), async, ϵ*t)
                end
            end
        end
    end
    tempnet.has_async_nodes = true
end

"""
Remove all asynchronus nodes of a TemporalGraph.

i.e. those put in place by add_async_nodes!
"""
function remove_async_nodes!(tempnet::QuNet.TemporalGraph)

    # No double dipping!
    if tempnet.has_async_nodes == false
        @warn("Redundent use of remove_async_nodes -- T.has_async_nodes == false
         Returning an unmodified TemporalGraph")
        return tempnet
    end

    offset = tempnet.nv * tempnet.steps
    for costkey in keys(zero_costvector())
        graph = tempnet.graph[costkey]
        # Assertion is no longer valid since not all async nodes from offset+1 to offset+N are generated
        # @assert nv(graph) == N * steps + N
        # Remove nodes with decreasing id so iteration works
        kill_list = reverse(collect((offset + 1) : (offset + tempnet.nv)))
        for i in kill_list
            rem_vertex!(graph, i)
        end
    end
    tempnet.has_async_nodes = false
end

"""
TODO: DOESN'T WORK

Remove the asynchronus edges of a TemporalGraph, excluding those connecting
specified layers
"""
# function remove_async_edges!(tempnet::QuNet.TemporalGraph, exclude_layers...)
#
#     @assert(all(0 < i <= tempnet.steps for i in exclude_layers))
#     # offset index for asynchronus nodes
#     offset = tempnet.nv * tempnet.steps
#
#     for costkey in keys(zero_costvector())
#         graph = tempnet.graph[costkey]
#         for async_node in offset + 1 : offset + tempnet.nv
#             nlist = all_neighbors(graph, async_node)
#             # remove nodes in nlist that are in the excluded layers
#             # BUG, this filters nodes, not edge layers.
#             filter!(x -> !((x-1) ÷ tempnet.nv in exclude_layers), nlist)
#             for node in nlist
#                 # Remove in both directions just in case
#                 rem_edge!(graph, async_node, node)
#                 rem_edge!(graph, node, async_node)
#             end
#         end
#     end
# end

"""
Convert temporal nodes into asynchronus nodes
If the node is already asynchronus, it is returned
Throws an exception if the node is out of bounds
"""
function temp_to_async(nv::Int, steps::Int, node::Int)
    offset = nv * steps
    if node <= 0 || node > offset + nv
        error("Node index out of bounds")
    elseif node > offset
        return node
    end
    mod = node % nv
    if mod == 0
        mod = nv
    end
    return mod + offset
end

"""
Get the temporal depth of a node from its index. If asynchronus, return -1
Throws an exception if the node is out of bounds.
"""
function node_timedepth(nv::Int, steps::Int, node::Int)
    offset = nv * steps
    if node <= 0 || node > offset + nv
        error("Node index out of bounds")
    elseif node > offset
        return -1
    end
    if node % nv == 0
         return (node-1) ÷ nv + 1
    end
    return (node ÷ nv) + 1
end

"""
Get the relative index of a node from its original index.
ie. for a graph with nv = 10, steps = 3, then node 14
is the same relative to node 4. Another way to think about this
is that it gets the index of the corresponding node in the top layer.
"""
function relative_index(nv::Int, steps::Int, node::Int)
    offset = nv * steps
    if node <= 0 || node > offset + nv
        error("Node index out of bounds")
    end
    mod = node % nv
    if mod == 0
        return nv
    end
    return mod
end

"""
For a given node at time T in a TemporalGraph, this function fixes its asynchronus
counterpart with respect to T This is done by removing all edges of the asynchronus node
except for the one connecting at time T.
"""
function fix_async_nodes_in_time!(tempnet::QuNet.TemporalGraph, node_list::Vector{Int})

    if tempnet.has_async_nodes == false
        return tempnet
    end

    offset = tempnet.steps * tempnet.nv
    # Check that nodes are valid non-temporal indices
    @assert(all(0 < i <= offset for i in node_list))

    for costkey in keys(zero_costvector())
        graph = tempnet.graph[costkey]
        for node in node_list
            # Relative_index of node
            rel = relative_index(tempnet.nv, tempnet.steps, node)
            # Get the timedepth of the node
            depth = node_timedepth(tempnet.nv, tempnet.steps, node)
            # Get the asynchronus counterpart of node
            async = temp_to_async(tempnet.nv, tempnet.steps, node)

            for i in 1:tempnet.steps
                if i != depth
                    hard_rem_edge!(graph, rel + (i-1)*tempnet.nv, async)
                end
            end
        end
    end
end

function gplot(temp_graph::TemporalGraph)
    gplot(temp_graph.graph["loss"], temp_graph.locs_x, temp_graph.locs_y, arrowlengthfrac=0.04)
end

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
export train!

# train core 1 {{{
function train_core1(p::Real,q::Real)
	return p * q
end

function train_core1{T<:Real}(p::Complex{T},q::Complex{T})
	return p * conj(q)
end

function train_core1{T<:Real}(p::Quaternion{T},q::Quaternion{T})
	return p * conj(q)
end
# }}}

# train core 2 {{{
function train_core2(x::Real,nn::Real)
	return x/nn
end

function train_core2{T<:Real}(z::Complex{T},nn::Real)
	return z/(2.0*nn)
end

function train_core2{T<:Real}(q::Quaternion{T},nn::Real)
	return q/(4.0*nn)
end
# }}}

function train!{T<:Number}(net::HopfieldNetwork, patterns::Array{T,2})
	nn,np = size(patterns)

	# calc synaptic weight (Hebbian learning rule)
	for p in 1:nn
		for q in p+1:nn
			s = zero(T)
			for mu in 1:np
				s += train_core1(patterns[p, mu], patterns[q, mu])
			end
			s = train_core2(s,nn)
			net.weight[p, q] = s
			net.weight[q, p] = conj(s)
		end
	end

	return
end

function train_core{T<:Real}(patterns::AbstractArray{T},nn::Real)
	tmp = BLAS.syrk('U','T',1.0/nn, patterns')
	for i=1:nn
		tmp[i,i] = zero(T)
		for j=i+1:nn
			tmp[j,i] = tmp[i,j]
		end
	end

	return tmp
end

function train_fast!{T <: Number}(net::HopfieldNetwork, patterns::Array{T})
	nn = length(net.state) # num of neurons
	np = size(patterns, 2) # num of patterns

	# calc synaptic weight
	net.Weight = train_core(patterns,nn)
	return
end


text:
struct ProductDistribution{D} <: AbstractMeasure
    marginals::Vector{D}
end

marginals(dist::ProductDistribution) = dist.marginals
marginal(dist::ProductDistribution, i::Integer) = dist.marginals[i]

length(dist::ProductDistribution) = length(marginals(dist))

## Simulation

function rand(rng::AbstractRNG, dist::ProductDistribution)
    return [rand(rng, marginal(dist, i)) for i = 1:length(dist)]
end

rand(dist::ProductDistribution) = rand(GLOBAL_RNG, dist)

## Likelihood

function logdensity(dist::ProductDistribution, x::AbstractVector)
    return sum(logdensity(marginal(dist, i), x[i]) for i = 1:length(dist))
end

## Prior likelihood

function logdensity(prior::ProductDistribution, dist::ProductDistribution)
    return sum(logdensity(marginal(prior, i), marginal(dist, i)) for i = 1:length(dist))
end

## Sufficient statistics

struct ProductSufficientStats{T<:Dists.SufficientStats}
    ss_marginals::Vector{T}
end

function suffstats(::Type{<:ProductDistribution{T}}, x::AbstractMatrix) where {T}
    return ProductSufficientStats([suffstats(T, x[i, :] for i in size(x, 1))])
end

function suffstats(
    ::Type{<:ProductDistribution{T}},
    x::AbstractMatrix,
    w::AbstractVector,
) where {T}
    return ProductSufficientStats([suffstats(T, x[i, :], w) for i in size(x, 1)])
end

function suffstats(
    ::Type{<:ProductDistribution{T}},
    prior::ProductDistribution{S},
    x::AbstractMatrix,
) where {T,S}
    return ProductSufficientStats([
        suffstats(T, marginal(prior, i), x[i, :]) for i in size(x, 1)
    ])
end

function suffstats(
    ::Type{<:ProductDistribution{T}},
    prior::ProductDistribution{S},
    x::AbstractMatrix,
    w::AbstractVector,
) where {T,S}
    return ProductSufficientStats([
        suffstats(T, marginal(prior, i), x[i, :], w) for i in size(x, 1)
    ])
end

## Fitting

function fit(pdtype::Type{<:ProductDistribution}, args...; kwargs...)
    return fit_mle(pdtype, args..., kwargs...)
end

function fit_mle(pdtype::Type{<:ProductDistribution}, args...; kwargs...)
    ss = suffstats(pdtype, args...; kwargs...)
    return fit_mle(pdtype, ss)
end


text:
using HypothesisTests, Base.Test
using HypothesisTests: default_tail

# data simulated under H_1
data_h1 = [
    -2.1844072302156263 1 -0.02187576598418861 1.4250842322231618 -2.0224314101114222
    -3.8237365225754556 1 -2.1844072302156263 -0.3986195023723653 -0.725318617093708
    -4.692622964314854 1 -3.8237365225754556 0.865359034508392 -0.9716762048860438
    -9.465165251695765 1 -4.692622964314854 -1.5248868209889381 -1.8178664695259354
    -3.784086411104036 1 -9.465165251695765 -0.018434833227298542 0.8719091772543835
    -7.682951955540945 1 -3.784086411104036 -0.6298048032223702 -1.6804917747830277
    -8.068822143998775 1 -7.682951955540945 0.3994394061296323 -0.5382489016869288
    -8.558608055778492 1 -8.068822143998775 -0.34662848289604303 -1.0915556882219897
    -10.905266364277344 1 -8.558608055778492 0.13174261149628444 -2.373156801968529
    -3.7593874757578343 1 -10.905266364277344 1.2873446803519781 1.3193356146665325
    -0.4071004208759922 1 -3.7593874757578343 0.6312912597626849 0.2928588585051935
    -0.21264428099357685 1 -0.4071004208759922 -1.2337261440146907 0.6677002797481758
    -0.2566792446386543 1 -0.21264428099357685 -0.8585845030636355 0.21112095190111557
    1.4652422469376796 1 -0.2566792446386543 0.2582660780400912 0.3254542082853489
    -0.03370549493934421 1 1.4652422469376796 -1.4599892950591262 0.3142665708225735
    -3.6655187094551787 1 -0.03370549493934421 0.5095274347206296 -1.2436129408018743
    -3.67007844594898 1 -3.6655187094551787 -1.2459007156231243 0.8724432818783148
    -4.269483828240617 1 -3.67007844594898 -1.3316534960103243 0.402333915295087
    -9.072321453339633 1 -4.269483828240617 -1.360732996599522 -0.918403552104298
    -6.298475056265913 1 -9.072321453339633 0.5429237557639659 0.564993224895932
    -8.851275959789886 1 -6.298475056265913 -2.63399370777994 -0.8745409882369648
    -7.240409805869097 1 -8.851275959789886 0.1499014406502077 -0.487019467609227
    -8.742193488440103 1 -7.240409805869097 -1.6587759684101517 -0.5287577864433531
    -8.085128207595446 1 -8.742193488440103 -2.748600587280482 0.6442930345934977
    -4.595734505262792 1 -8.085128207595446 -0.4180662233306449 0.909924517665989
    -8.876081007098774 1 -4.595734505262792 -1.4313262030289144 -1.3628300743751958
    -6.097510119608872 1 -8.876081007098774 1.1540597669111703 0.8739856083349014
    -8.187395934761742 1 -6.097510119608872 -0.4855284875097715 -0.8548933573000406
    -4.630296021087888 1 -8.187395934761742 0.7754380459705584 0.7550975795526141
    -6.678197570842665 1 -4.630296021087888 0.7560596380052317 -0.8436185092902743
    -8.454654098952108 1 -6.678197570842665 -1.626729284162485 0.27766300810348893
    -8.87315711823648 1 -8.454654098952108 -0.944684415100926 -0.9030356944859604
    -8.299137304666294 1 -8.87315711823648 0.658322671059814 -1.513350928595663
    -7.66583122908761 1 -8.299137304666294 -0.4489512342884181 -0.6593759619410697
    -5.544754187110598 1 -7.66583122908761 1.1851174294240252 -0.7403475629601508
    -0.2346062802414116 1 -5.544754187110598 1.9716204287413468 1.1351529072038815
    -2.7080211703552033 1 -0.2346062802414116 -0.0740194597978335 -0.8957439911719416
    -4.405304542289478 1 -2.7080211703552033 -0.20212227601885754 0.02136996876872505
    -10.874421844205603 1 -4.405304542289478 -1.458002434697322 -1.5544394018554584
    -7.980300798148038 1 -10.874421844205603 1.797342002330958 0.193735930526397
    -6.092985977023801 1 -7.980300798148038 -0.9064866353733989 0.9206450790295183
    -5.429540806278794 1 -6.092985977023801 -0.9702203522809195 0.16364495372354845
    -2.8135611982312048 1 -5.429540806278794 -0.7689287953327062 0.2606593892993496
    -5.233993628762067 1 -2.8135611982312048 -0.3959933139142994 -2.620564252102542
    -4.636500469478656 1 -5.233993628762067 0.4751826349484905 -1.0562153109753674
    -3.951599322839157 1 -4.636500469478656 0.6386323305899096 -0.0903512992302431
    -4.5338888862750215 1 -3.951599322839157 -0.3987477530370317 0.27865272355818926
    -3.27448245549044 1 -4.5338888862750215 -0.0803466446475703 -0.06596847546469198
    -2.116669048507199 1 -3.27448245549044 1.0018399414103412 -0.1569479584427068
    -5.632636058340474 1 -2.116669048507199 -0.7048446599977012 -1.8390689017857673
    -6.612580447433536 1 -5.632636058340474 -1.4340505462257125 -0.5674519382746356
    -6.047621290327909 1 -6.612580447433536 0.4347509475695458 -0.9192186367107573
    -12.316195547709809 1 -6.047621290327909 -0.09220577556331351 -3.8819978594268876
    -9.790653346165218 1 -12.316195547709809 -0.46408289171823325 0.8665458026221579
    -7.081636830140645 1 -9.790653346165218 2.187832361106393 0.432306709021104
    -7.798466676804453 1 -7.081636830140645 0.5791414166174949 -0.3229959302497174
    -8.48368570601047 1 -7.798466676804453 -0.9147242711971281 0.3402161710448687
    -6.844741808358232 1 -8.48368570601047 -0.1550598509968508 -0.29961787965293785
    -1.9143881179353244 1 -6.844741808358232 0.060717820312243403 1.9136035721789713
    -6.350825197915615 1 -1.9143881179353244 -0.27085887966353506 -1.6190985435203415
    -5.386740290175477 1 -6.350825197915615 0.11104289827811374 0.5169712744781768
    -8.201462522665064 1 -5.386740290175477 0.6427956685597893 -1.2825811121125588
    -6.86037895816535 1 -8.201462522665064 -1.0932960849422242 1.0840840119878734
    -7.681177025769711 1 -6.86037895816535 0.0867583435648074 -0.7944197379913668
    -8.700601165824938 1 -7.681177025769711 -0.07077890859801883 -0.025668064894330093
    -6.737037436370262 1 -8.700601165824938 1.2497099131508578 -0.011462506736218826
    -3.6689503491552253 1 -6.737037436370262 0.1512033860208509 1.4678868411216563
    -1.249607364746875 1 -3.6689503491552253 0.35515785025036173 1.315845372889509
    -1.9479972265280567 1 -1.249607364746875 -0.9253434974479612 0.27236940715678176
    0.06459657520168782 1 -1.9479972265280567 -0.233103327838231 0.44237536440754566
    1.6736047640968876 1 0.06459657520168782 1.3606109850208474 -1.1966483101519632
    3.3481475962652025 1 1.6736047640968876 0.7311973736313007 0.48457993854995624
    4.395242135786573 1 3.3481475962652025 0.4787377009574211 1.4147590837932418
    5.718254592907931 1 4.395242135786573 0.6051343777576946 1.0451679003902126
    7.701764680294246 1 5.718254592907931 1.276972611856289 -0.1978660512585745
    8.761139787455583 1 7.701764680294246 -1.3961992590989196 0.32892107972067225
    8.130728492876539 1 8.761139787455583 0.11709983671299977 -1.7983067508193347
    12.0528015874146 1 8.130728492876539 0.6266845348314155 1.44050014454746
    11.188183563801353 1 12.0528015874146 0.3133823643278484 -0.5232908450257276
    12.92107497611711 1 11.188183563801353 0.9862522782876494 0.43490251042216804
    12.271768495692909 1 12.92107497611711 0.5066085599761962 0.07774415899993356
    14.189399733550822 1 12.271768495692909 0.8718155606585614 0.14324697267944667
    15.088143644678663 1 14.189399733550822 0.6995645399732705 -0.43422769970595093
    17.365142073656695 1 15.088143644678663 2.0096551932105817 0.24387704638876334
    15.474852839205559 1 17.365142073656695 0.15269868152330113 -0.12389276233418343
    14.228935730675236 1 15.474852839205559 -0.7879019375423509 0.42889110309176043
    17.394434630323023 1 14.228935730675236 0.9986772266084378 1.189248289168407
    16.54384322429009 1 17.394434630323023 -0.4953901743665539 0.2059465326047665
    18.353069936157006 1 16.54384322429009 0.3638388194911837 1.0218854109275881
    16.07119383937772 1 18.353069936157006 0.6753443406355725 -1.1990648633034413
    20.977942062323052 1 16.07119383937772 0.6008899714463747 2.3805130457880845
    23.626966845968067 1 20.977942062323052 -1.5021368491168559 2.1896086828053907
    27.902989083231745 1 23.626966845968067 -0.3727356407152467 2.1299921566412503
    24.021974044852822 1 27.902989083231745 -1.2766407160812419 -0.3876601960862665
    22.470178938207653 1 24.021974044852822 0.6451054034379425 -0.23211025519653547
    22.71778230013325 1 22.470178938207653 1.9058844942684128 0.20626950528981075
    18.975581322758284 1 22.71778230013325 -1.5881288112102452 1.107610369079468
    19.971416875068105 1 18.975581322758284 0.29738396177383514 0.8667513992731141
    23.11402529567158 1 19.971416875068105 0.3028824231221631 1.6776251628069427
    17.780386660161668 1 23.11402529567158 -0.10176910274960058 -1.4298441320561919
]

coeff_vec = data_h1[:,2:end]\data_h1[:,1]
res_vec = data_h1[:,1] - data_h1[:,2:end]*coeff_vec

t = BreuschGodfreyTest(data_h1[:,2:end],res_vec,4)

@test t.n == 100
@test t.lag == 4
@test t.BG ≈ 31.39810637185552
@test pvalue(t) ≈ 2.5390992557054064e-6
@test default_tail(t) == :right
show(IOBuffer(), t)

t = BreuschGodfreyTest(data_h1[:,2:end],res_vec,2,false)

@test t.n == 98
@test t.lag == 2
@test t.BG ≈ 27.71445552148237
@test pvalue(t) ≈ 9.591409846298036e-7
show(IOBuffer(), t)

coeff_vec = data_h1[:,3:end]\data_h1[:,1] # no constant
res_vec = data_h1[:,1] - data_h1[:,3:end]*coeff_vec

t = BreuschGodfreyTest(data_h1[:,3:end],res_vec,10)

@test t.n == 100
@test t.lag == 10
@test t.BG ≈ 34.562955523063515
@test pvalue(t) ≈ 0.0001482189344859451
show(IOBuffer(), t)

# data simulated under H_0
data_h0 = [
    -2.6393758379263925 1 -0.22997735515901174 1.4250842322231618 -2.0224314101114222
    -3.8148641358053834 1 -2.6393758379263925 -0.3986195023723653 -0.725318617093708
    -4.677326144975626 1 -3.8148641358053834 0.865359034508392 -0.9716762048860438
    -8.864079316218326 1 -4.677326144975626 -1.5248868209889381 -1.8178664695259354
    -3.460887006123175 1 -8.864079316218326 -0.018434833227298542 0.8719091772543835
    -8.998633967744288 1 -3.460887006123175 -0.6298048032223702 -1.6804917747830277
    -7.665975510821412 1 -8.998633967744288 0.3994394061296323 -0.5382489016869288
    -8.56702198923559 1 -7.665975510821412 -0.34662848289604303 -1.0915556882219897
    -11.214833972884769 1 -8.56702198923559 0.13174261149628444 -2.373156801968529
    -4.446594124802662 1 -11.214833972884769 1.2873446803519781 1.3193356146665325
    -1.6998765285430706 1 -4.446594124802662 0.6312912597626849 0.2928588585051935
    -1.6972964464315492 1 -1.6998765285430706 -1.2337261440146907 0.6677002797481758
    -1.2272304633377735 1 -1.6972964464315492 -0.8585845030636355 0.21112095190111557
    0.30586906328114577 1 -1.2272304633377735 0.2582660780400912 0.3254542082853489
    -1.041825499276659 1 0.30586906328114577 -1.4599892950591262 0.3142665708225735
    -4.0872691626960425 1 -1.041825499276659 0.5095274347206296 -1.2436129408018743
    -2.8079705764801184 1 -4.0872691626960425 -1.2459007156231243 0.8724432818783148
    -3.5781229211748493 1 -2.8079705764801184 -1.3316534960103243 0.402333915295087
    -7.919753510576073 1 -3.5781229211748493 -1.360732996599522 -0.918403552104298
    -4.050334580805227 1 -7.919753510576073 0.5429237557639659 0.564993224895932
    -7.385113567278499 1 -4.050334580805227 -2.63399370777994 -0.8745409882369648
    -6.389907039488093 1 -7.385113567278499 0.1499014406502077 -0.487019467609227
    -8.362337662947732 1 -6.389907039488093 -1.6587759684101517 -0.5287577864433531
    -7.698895749951619 1 -8.362337662947732 -2.748600587280482 0.6442930345934977
    -5.183985240958993 1 -7.698895749951619 -0.4180662233306449 0.909924517665989
    -9.354302625888938 1 -5.183985240958993 -1.4313262030289144 -1.3628300743751958
    -6.0060845548714905 1 -9.354302625888938 1.1540597669111703 0.8739856083349014
    -7.251074336523842 1 -6.0060845548714905 -0.4855284875097715 -0.8548933573000406
    -3.750382446404393 1 -7.251074336523842 0.7754380459705584 0.7550975795526141
    -5.770679903165536 1 -3.750382446404393 0.7560596380052317 -0.8436185092902743
    -6.1571290241435745 1 -5.770679903165536 -1.626729284162485 0.27766300810348893
    -6.758646960354187 1 -6.1571290241435745 -0.944684415100926 -0.9030356944859604
    -7.251801456499875 1 -6.758646960354187 0.658322671059814 -1.513350928595663
    -7.120615787554741 1 -7.251801456499875 -0.4489512342884181 -0.6593759619410697
    -5.588509539557388 1 -7.120615787554741 1.1851174294240252 -0.7403475629601508
    -0.5390053133408673 1 -5.588509539557388 1.9716204287413468 1.1351529072038815
    -2.5300623604990955 1 -0.5390053133408673 -0.0740194597978335 -0.8957439911719416
    -3.7798892423452592 1 -2.5300623604990955 -0.20212227601885754 0.02136996876872505
    -9.04942964730818 1 -3.7798892423452592 -1.458002434697322 -1.5544394018554584
    -5.4354744982455365 1 -9.04942964730818 1.797342002330958 0.193735930526397
    -3.5679532514488668 1 -5.4354744982455365 -0.9064866353733989 0.9206450790295183
    -3.2063386742141633 1 -3.5679532514488668 -0.9702203522809195 0.16364495372354845
    -1.051934916538913 1 -3.2063386742141633 -0.7689287953327062 0.2606593892993496
    -4.863145941477073 1 -1.051934916538913 -0.3959933139142994 -2.620564252102542
    -5.212678725939295 1 -4.863145941477073 0.4751826349484905 -1.0562153109753674
    -4.614614610635563 1 -5.212678725939295 0.6386323305899096 -0.0903512992302431
    -4.424400676415464 1 -4.614614610635563 -0.3987477530370317 0.27865272355818926
    -2.5809040973960595 1 -4.424400676415464 -0.0803466446475703 -0.06596847546469198
    -2.283035131341285 1 -2.5809040973960595 1.0018399414103412 -0.1569479584427068
    -4.838392113346874 1 -2.283035131341285 -0.7048446599977012 -1.8390689017857673
    -6.765964527239637 1 -4.838392113346874 -1.4340505462257125 -0.5674519382746356
    -6.146713996908273 1 -6.765964527239637 0.4347509475695458 -0.9192186367107573
    -12.93695213339396 1 -6.146713996908273 -0.09220577556331351 -3.8819978594268876
    -10.328147452628222 1 -12.93695213339396 -0.46408289171823325 0.8665458026221579
    -7.312623658528083 1 -10.328147452628222 2.187832361106393 0.432306709021104
    -6.601242580580354 1 -7.312623658528083 0.5791414166174949 -0.3229959302497174
    -7.007982630266133 1 -6.601242580580354 -0.9147242711971281 0.3402161710448687
    -4.711869773110432 1 -7.007982630266133 -0.1550598509968508 -0.29961787965293785
    -1.3112388084541342 1 -4.711869773110432 0.060717820312243403 1.9136035721789713
    -4.778553640411378 1 -1.3112388084541342 -0.27085887966353506 -1.6190985435203415
    -3.597048209305651 1 -4.778553640411378 0.11104289827811374 0.5169712744781768
    -5.91625887800437 1 -3.597048209305651 0.6427956685597893 -1.2825811121125588
    -3.8341200305975844 1 -5.91625887800437 -1.0932960849422242 1.0840840119878734
    -5.05155415725161 1 -3.8341200305975844 0.0867583435648074 -0.7944197379913668
    -5.902658917952474 1 -5.05155415725161 -0.07077890859801883 -0.025668064894330093
    -3.0148968638612246 1 -5.902658917952474 1.2497099131508578 -0.011462506736218826
    -0.5435623298882797 1 -3.0148968638612246 0.1512033860208509 1.4678868411216563
    2.576402881110434 1 -0.5435623298882797 0.35515785025036173 1.315845372889509
    1.840001826392218 1 2.576402881110434 -0.9253434974479612 0.27236940715678176
    3.7629568916181633 1 1.840001826392218 -0.233103327838231 0.44237536440754566
    4.300779090152549 1 3.7629568916181633 1.3606109850208474 -1.1966483101519632
    4.890327561532721 1 4.300779090152549 0.7311973736313007 0.48457993854995624
    6.712083346070772 1 4.890327561532721 0.4787377009574211 1.4147590837932418
    8.918662089794555 1 6.712083346070772 0.6051343777576946 1.0451679003902126
    10.889228410523366 1 8.918662089794555 1.276972611856289 -0.1978660512585745
    10.771798896063286 1 10.889228410523366 -1.3961992590989196 0.32892107972067225
    8.897958495974741 1 10.771798896063286 0.11709983671299977 -1.7983067508193347
    11.231705359232382 1 8.897958495974741 0.6266845348314155 1.44050014454746
    11.085832271063609 1 11.231705359232382 0.3133823643278484 -0.5232908450257276
    12.02548460222439 1 11.085832271063609 0.9862522782876494 0.43490251042216804
    11.819287696755126 1 12.02548460222439 0.5066085599761962 0.07774415899993356
    13.96911393915951 1 11.819287696755126 0.8718155606585614 0.14324697267944667
    13.840528753974972 1 13.96911393915951 0.6995645399732705 -0.43422769970595093
    15.684175759900603 1 13.840528753974972 2.0096551932105817 0.24387704638876334
    14.08213146714004 1 15.684175759900603 0.15269868152330113 -0.12389276233418343
    13.303842222767337 1 14.08213146714004 -0.7879019375423509 0.42889110309176043
    16.40951365835989 1 13.303842222767337 0.9986772266084378 1.189248289168407
    15.406011864201629 1 16.40951365835989 -0.4953901743665539 0.2059465326047665
    17.105004411081747 1 15.406011864201629 0.3638388194911837 1.0218854109275881
    14.766483437171006 1 17.105004411081747 0.6753443406355725 -1.1990648633034413
    19.48191551080787 1 14.766483437171006 0.6008899714463747 2.3805130457880845
    22.133552497574904 1 19.48191551080787 -1.5021368491168559 2.1896086828053907
    25.492693395853813 1 22.133552497574904 -0.3727356407152467 2.1299921566412503
    20.970851703684147 1 25.492693395853813 -1.2766407160812419 -0.3876601960862665
    19.83986700193583 1 20.970851703684147 0.6451054034379425 -0.23211025519653547
    20.241165411067264 1 19.83986700193583 1.9058844942684128 0.20626950528981075
    17.31664682866259 1 20.241165411067264 -1.5881288112102452 1.107610369079468
    19.675283574585194 1 17.31664682866259 0.29738396177383514 0.8667513992731141
    21.79754881195714 1 19.675283574585194 0.3028824231221631 1.6776251628069427
    16.68579877178968 1 21.79754881195714 -0.10176910274960058 -1.4298441320561919
]

coeff_vec = data_h0[:,2:end]\data_h0[:,1]
res_vec = data_h0[:,1] - data_h0[:,2:end]*coeff_vec

t = BreuschGodfreyTest(data_h0[:,2:end],res_vec,4)

@test t.n == 100
@test t.lag == 4
@test t.BG ≈ 1.9781996740849623
@test pvalue(t) ≈ 0.7397687483602012
show(IOBuffer(), t)

t = BreuschGodfreyTest(data_h0[:,2:end],res_vec,2,false)

@test t.n == 98
@test t.lag == 2
@test t.BG ≈ 0.03720525092311448
@test pvalue(t) ≈ 0.9815693354165801
show(IOBuffer(), t)

coeff_vec = data_h0[:,3:end]\data_h0[:,1] # no constant
res_vec = data_h0[:,1] - data_h0[:,3:end]*coeff_vec

t = BreuschGodfreyTest(data_h0[:,3:end],res_vec,10)

@test t.n == 100
@test t.lag == 10
@test t.BG ≈ 11.156984521656943
@test pvalue(t) ≈ 0.3454190353708124
show(IOBuffer(), t)


text:
module GridapTests

using Test

@time @testset "Poisson" begin include("PoissonTests.jl") end

@time @testset "PoissonDG" begin include("PoissonDGTests.jl") end

@time @testset "Poisson1D" begin include("Poisson1DTests.jl") end

@time @testset "PLaplacian" begin include("PLaplacianTests.jl") end

@time @testset "PLaplacianWithAutodiff" begin include("PLaplacianWithAutodiffTests.jl") end

@time @testset "StokesTaylorHood" begin include("StokesTaylorHoodTests.jl") end

@time @testset "StokesDG" begin include("StokesDGTests.jl") end

@time @testset "StokesNitsche" begin include("StokesNitscheTests.jl") end

@time @testset "Darcy" begin include("DarcyTests.jl") end

@time @testset "PeriodicDarcy" begin include("PeriodicDarcyTests.jl") end

@time @testset "SurfaceCoupling" begin include("SurfaceCouplingTests.jl") end

@time @testset "IsotropicDamage" begin include("IsotropicDamageTests.jl") end

@time @testset "PhysicalPoisson" begin include("PhysicalPoissonTests.jl") end

end # module


text:
### Composite Curve ###

mutable struct CompositeIRCurve <: AbstractIRCurve
    date::Date
    list::Vector{AbstractIRCurve}

    function CompositeIRCurve(curve_tuple::AbstractIRCurve...)
        @assert !isempty(curve_tuple) "Empty list of curves for CompositeIRCurve"
        date = curve_get_date(curve_tuple[1])
        len = length(curve_tuple)
        if len > 1
            for i in 2:len
                @assert date == curve_get_date(curve_tuple[i]) "All curve dates should match for CompositeIRCurve"
            end
        end

        list = Vector{AbstractIRCurve}(undef, len)
        for i in 1:len
            @inbounds list[i] = curve_tuple[i]
        end

        new(date, list)
    end
end

curve_get_date(c::CompositeIRCurve) = c.date

function ERF(curve::CompositeIRCurve, maturity::Date)
    local _erf::Float64 = 1.0
    for c in curve.list
        _erf *= ERF(c, maturity)
    end
    return _erf
end

# Unoptimized vector function for ERF
function ERF(curve::CompositeIRCurve, maturity_vec::Vector{Date})
    len = length(maturity_vec)
    result = Vector{Float64}(undef, len)
    for i in 1:len
        @inbounds result[i] = ERF(curve, maturity_vec[i])
    end
    return result
end

# Unoptimized vector function for discountfactor
function discountfactor(curve::CompositeIRCurve, maturity_vec::Vector{Date})
    len = length(maturity_vec)
    result = Vector{Float64}(undef, len)
    for i in 1:len
        @inbounds result[i] = discountfactor(curve, maturity_vec[i])
    end
    return result
end

ERF_to_rate(curve::CompositeIRCurve, ERF::Float64, t::Float64) = error("Function not available for CompositeIRCurve")
forward_rate(curve::CompositeIRCurve, forward_date::Date, maturity::Date) = error("Function not available for CompositeIRCurve")
zero_rate(curve::CompositeIRCurve, maturity::Date) = error("Function not available for CompositeIRCurve")
zero_rate(curve::CompositeIRCurve, maturity_vec::Vector{Date}) = error("Function not available for CompositeIRCurve")


text:
@testsuite "constructors" (AT, eltypes)->begin
    @testset "direct" begin
        for T in eltypes
            B = AT{T}(undef, 10)
            @test B isa AT{T,1}
            @test size(B) == (10,)
            @test eltype(B) == T

            B = AT{T}(undef, 10, 10)
            @test B isa AT{T,2}
            @test size(B) == (10, 10)
            @test eltype(B) == T

            B = AT{T}(undef, (10, 10))
            @test B isa AT{T,2}
            @test size(B) == (10, 10)
            @test eltype(B) == T
        end

        # compare against Array
        for typs in [(), (Int,), (Int,1), (Int,2), (Float32,), (Float32,1), (Float32,2)],
            args in [(), (1,), (1,2), ((1,),), ((1,2),),
                    (undef,), (undef, 1,), (undef, 1,2), (undef, (1,),), (undef, (1,2),),
                    (Int,), (Int, 1,), (Int, 1,2), (Int, (1,),), (Int, (1,2),),
                    ([1,2],), ([1 2],)]
            cpu = try
                Array{typs...}(args...)
            catch ex
                isa(ex, MethodError) || rethrow()
                nothing
            end

            gpu = try
                AT{typs...}(args...)
            catch ex
                isa(ex, MethodError) || rethrow()
                cpu == nothing || rethrow()
                nothing
            end

            if cpu == nothing
                @test gpu == nothing
            else
                @test typeof(cpu) == typeof(convert(Array, gpu))
            end
        end
    end

    @testset "similar" begin
        for T in eltypes
            B = AT{T}(undef, 10)

            B = similar(B, Int32, 11, 15)
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(B, T)
            @test B isa AT{T,2}
            @test size(B) == (11, 15)
            @test eltype(B) == T

            B = similar(B, (5,))
            @test B isa AT{T,1}
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(B, 7)
            @test B isa AT{T,1}
            @test size(B) == (7,)
            @test eltype(B) == T

            B = similar(AT{Int32}, (11, 15))
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(AT{T}, (5,))
            @test B isa AT{T,1}
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(AT{T}, 7)
            @test B isa AT{T,1}
            @test size(B) == (7,)
            @test eltype(B) == T

            B = similar(Broadcast.Broadcasted(*, (B, B)), T)
            @test B isa AT{T,1}
            @test size(B) == (7,)
            @test eltype(B) == T

            B = similar(Broadcast.Broadcasted(*, (B, B)), Int32, (11, 15))
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32
        end
    end

    @testset "convenience" begin
        for T in eltypes
            A = AT(rand(T, 3))
            b = rand(T)
            fill!(A, b)
            @test A isa AT{T,1}
            @test Array(A) == fill(b, 3)

            A = zero(AT(rand(T, 2)))
            @test A isa AT{T,1}
            @test Array(A) == zero(rand(T, 2))

            A = zero(AT(rand(T, 2, 2)))
            @test A isa AT{T,2}
            @test Array(A) == zero(rand(T, 2, 2))

            A = zero(AT(rand(T, 2, 2, 2)))
            @test A isa AT{T,3}
            @test Array(A) == zero(rand(T, 2, 2, 2))

            A = one(AT(rand(T, 2, 2)))
            @test A isa AT{T,2}
            @test Array(A) == one(rand(T, 2, 2))

            A = oneunit(AT(rand(T, 2, 2)))
            @test A isa AT{T,2}
            @test Array(A) == oneunit(rand(T, 2, 2))
        end
    end

    @testset "conversions" begin
        for T in eltypes
            Bc = round.(rand(10, 10) .* 10.0)
            B = AT{T}(Bc)
            @test size(B) == (10, 10)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(T, 10)
            B = AT(Bc)
            @test size(B) == (10,)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(T, 10, 10)
            B = AT{T, 2}(Bc)
            @test size(B) == (10, 10)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            intervals = Dict(
                Float16 => -2^11:2^11,
                Float32 => -2^24:2^24,
                Float64 => -2^53:2^53,
            )

            Bc = rand(Int8, 3, 3, 3)
            B = convert(AT{T, 3}, Bc)
            @test size(B) == (3, 3, 3)
            @test eltype(B) == T
            @test Array(B) ≈ Bc
        end
    end

    @testset "uniformscaling" begin
        for T in eltypes
            x = Matrix{T}(I, 4, 2)

            x1 = AT{T, 2}(I, 4, 2)
            x2 = AT{T}(I, (4, 2))
            x3 = AT{T, 2}(I, (4, 2))

            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x

            x = Matrix(T(3) * I, 2, 4)
            x1 = AT(T(3) * I, 2, 4)
            @test eltype(x1) == T
            @test Array(x1) ≈ x
        end
    end
end


text:
"""
    mblevitus(cmd0::String=""; kwargs...)

Create a water velocity profile which is representative of the mean annual water column for a specified
1 degree by 1 degree region.

Parameters
----------

- **A** | **all4** :: [Type => Bool]

    Pint also depth, velocity, temperature, salinity.
- **L** | **location** :: [Type => Str | Tuple]		``Arg = lon/lat``

    Sets the longitude and latitude of the location of the water velocity profile.
- **O** | **outfile** | **out_file** :: [Type => Str]

    Write the SVP to <outfile>.
- **H** | **help** :: [Type => Bool]

    Print out program's description.
- **z** | **z_down** :: [Type => Bool]

    Makes Z axes positive down (default here is Z-up).
- $(GMT.opt_V)
"""
function mblevitus(cmd0::String=""; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("mblevitus", cmd0)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd = parse_common_opts(d, "", [:yx :V_params :o])[1]
	cmd = add_opt(d, cmd, "L", [:L :location :R])
	(!occursin("-L", cmd)) && (cmd *= " -L0/0")
	cmd = parse_these_opts(cmd, d, [[:A :all4], [:O :outfile :output_file], [:H :help], [:z :z_down]])
	((val = find_in_dict(d, [:z :z_down])[1]) === nothing) && (cmd *= " -z")	# Means here default is Z-up
	(!occursin("-:", cmd)) && (cmd = "-:" * cmd)		# Means that the default here (contrary to C version) is speed-dept

	finish_PS_module(d, "mblevitus " * cmd, "", true, false, false)
end


text:
using DelimitedFiles, Unicode

"""

# read_variational

Read variational sample output files created by cmdstan.

### Method
```julia
read_variational(m::Stanmodel)
```

### Required arguments
```julia
* `m::Stanmodel`    : Stanmodel object
```

"""
function read_variational(m::Stanmodel)

  local a3d, monitors, index, idx, indvec, ftype

  ftype = lowercase(string(typeof(m.method)))

  for i in 1:m.nchains
    if isfile("$(m.name)_$(ftype)_$(i).csv")
      open("$(m.name)_$(ftype)_$(i).csv") do instream
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)
        idx = split(strip(line), ",")
        index = [idx[k] for k in 1:length(idx)]
        if length(m.monitors) == 0
          indvec = 1:length(index)
        else
          indvec = findall((in)(m.monitors), index)
        end
        if i == 1
          a3d = fill(0.0, m.method.output_samples, length(indvec), m.nchains)
        end
        skipchars(isspace, instream, linecomment='#')
        for j in 1:m.method.output_samples
          skipchars(isspace, instream, linecomment='#')
          line = Unicode.normalize(readline(instream), newline2lf=true)
          if eof(instream) && length(line) < 2
            return
          else
            flds = parse.(Float64, (split(strip(line), ",")))
            flds = reshape(flds[indvec], 1, length(indvec))
            a3d[j,:,i] = flds
          end
        end
      end
    end
  end

  cnames = convert.(String, idx[indvec])

  (a3d, cnames)

end


text:
# Helper routines

real(x...) = Base.real(x...)
real{T<:Real}(::Type{T}) = T
real{T<:Real}(::Type{Complex{T}}) = T

eps{T<:Real}(::Type{T}) = Base.eps(T)
eps{T<:Real}(::Type{Complex{T}}) = eps(real(T))

#
# This function provides a convenient way to query or specify the BigFloat precision.
#
digits(n::Int) = setprecision(round(Int,ceil(n*log2(10))))
digits() = round(Int,floor(precision(BigFloat)*log10(2)))

h{T<:Number}(::Type{T},n::Int) = log(convert(T,π)*n)/n
h(n::Int) = h(Float64,n)

ω{T<:Number}(β::T,t::T) = exp(-β*(cosh(t)-one(T)))

function interlace{T<:Number}(u::Vector{T})
    n = length(u)
    v = zeros(T,n)
    if isodd(n)
        m = div(n,2) # n = 2m+1
        v[1] = u[1]
        for i = 1:m
          v[2i],v[2i+1] = u[i+m+1],u[i+1]
        end
    else
        m = div(n,2) # n = 2m
        for i=1:m
            v[2i-1],v[2i] = u[i+m],u[i]
        end
    end
    v
end

function interlace2{T<:Number}(u::Vector{T})
    n = length(u)
    v = zeros(T,n)
    if isodd(n)
        m = div(n,2) # n = 2m+1
        v[1] = u[m+1]
        for i = 1:m
          v[2i],v[2i+1] = u[m+i+1],u[m-i+1]
        end
    else
        m = div(n,2) # n = 2m
        for i=1:m
            v[2i-1],v[2i] = u[m-i+1],u[m+i]
        end
    end
    v
end


text:
"""
================================================================================
# Generic support types and functions
================================================================================
"""

export read_file, get_json, PKG_ROOT_DIR, check_attributes, generate, join_str, set_random_id, show_dict, funcName, getFunc, get_attributes

# read string from file
function read_file(name::String)
  content = string()
  if isfile(name) == true
    content = open(name) do file
     read(file, String)
    end
  else content = false
  end
  return content
end

 # Get dict from json file
 function get_json(file::String)
   read_file(file) != false ? JSON.parse(read_file(file)) : false
 end

 # Check for mandatory attributes
 function check_attributes(attributes::Dict,default_attributes::Dict)
       for key in keys(default_attributes)
           attributes[key] = haskey(attributes,key) ? attributes[key] : default_attributes[key]
       end
       return attributes
 end

 # Generate functions for language-specific elements
 function generate(lib::Dict,pref::String=string())
   for (e,v) in lib
       # Prefix is optional, so if no prefix is given, just use the element name
       name = !isempty(pref) ? string(pref,e) : e
       fn = Symbol(name)
       default = !isempty(v) ? v : Dict()
       eval(
         quote
              $fn = function (attributes::Dict=$default,children::String=string())
                     # attributes = check_attributes(attributes,$default)
                     return !isempty(children) ? element($e,attributes,children) : element($e,attributes)
                   end
              export $fn
         end
       )
   end
 end

 # Helper function to join an array of strings into one string, i.e. viewBox attribute
 function join_str(args)
   chop(join(map(v -> string("$v "), args)),head = 0, tail = 1)
 end

 # Set random id values for elements.
 # Alphanumeric characters only because href values with numbers as the first character do not work in svg documents
 function set_random_id(l::Int64=8)
   Random.randstring('a':'z', l)
 end

 # Show key and values of a dictionary as a string list
 function show_dict(d::Dict)
   content = string()
   for (key,value) in sort(collect(d))
         content = (string(content, "- ", key," : ",value, "\n"))
     end
     return """
            $content
            """
 end

 # Return function name of the current function
 # Reference: https://discourse.julialang.org/t/how-to-print-function-name-and-source-file-line-number/43486
 macro funcName()
    quote
         st = stacktrace(backtrace())
         myf = ""
         for frm in st
             funcname = frm.func
             if frm.func != :backtrace && frm.func!= Symbol("macro expansion")
                 myf = frm.func
                 break
             end
         end
         string(myf)
     end
 end

 # turn string into functions
 function getFunc(name)
     fn = Symbol(name)
     eval(
        quote
            return $fn()
        end
     )
 end

 # Return a copy of a dict value
 function get_attributes(d::Dict, value::String)
   return(copy(d[value]))
 end


text:
function bernoulli(n)
    A = Vector{Rational{BigInt}}(undef, n + 1)
    for m = 0 : n
        A[m + 1] = 1 // (m + 1)
        for j = m : -1 : 1
            A[j] = j * (A[j] - A[j + 1])
        end
    end
    return A[1]
end

function display(n)
    B = map(bernoulli, 0 : n)
    pad = mapreduce(x -> ndigits(numerator(x)) + Int(x < 0), max, B)
    argdigits = ndigits(n)
    for i = 0 : n
        if numerator(B[i + 1]) & 1 == 1
            println(
                "B(", lpad(i, argdigits), ") = ",
                lpad(numerator(B[i + 1]), pad), " / ", denominator(B[i + 1])
            )
        end
    end
end

display(60)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ---
# title: 1026. Maximum Difference Between Node and Ancestor
# id: problem1026
# author: Indigo
# date: 2022-04-14
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/maximum-difference-between-node-and-ancestor/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, find the maximum value `V` for which there
# exist **different** nodes `A` and `B` where `V = |A.val - B.val|` and `A` is
# an ancestor of `B`.
#
# A node `A` is an ancestor of `B` if either: any child of `A` is equal to `B`,
# or any child of `A` is an ancestor of `B`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/09/tmp-tree.jpg)
#
#
#
#     Input: root = [8,3,10,1,6,null,14,null,null,4,7,13]
#     Output: 7
#     Explanation: We have various ancestor-node differences, some of which are given below :
#     |8 - 3| = 5
#     |3 - 7| = 4
#     |8 - 1| = 7
#     |10 - 13| = 3
#     Among all possible differences, the maximum value of 7 is obtained by |8 - 1| = 7.
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/09/tmp-tree-1.jpg)
#
#
#
#     Input: root = [1,null,2,null,0,3]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[2, 5000]`.
#   * `0 <= Node.val <= 105`
#
#
## @lc code=start
using LeetCode

function max_ancestor_diff(root::TreeNode{Int})
    res = 0
    function min_max_descendant(root::TreeNode{Int})
        minres, maxres = root.val, root.val
        for child in (:left, :right)
            cd_node = getproperty(root, child)
            if !isnothing(cd_node)
                minl, maxl = min_max_descendant(cd_node)
                res = max(res, abs(root.val - minl), abs(root.val - maxl))
                minres = min(minres, minl)
                maxres = max(maxres, maxl)
            end
        end
        return minres, maxres
    end
    min_max_descendant(root)
    return res
end
## @lc code=end


text:
test_path = joinpath(Pkg.dir("SixelGraphics"),"test")
for test_file in readdir(test_path)
    include(joinpath(test_path, test_file))
end


text:
#= 本プログラムファイルがあるフォルダ内にある、TLEを入力したtxtファイルを読み込み、
   TLEの要素ごとに分割するプログラム                                           =#

txtName="ISS_TLE.txt"
open("ISS_TLE.txt","r") do fp
    lines = readlines(fp)
    line1 = split(lines[2])
    line2 = split(lines[3])

    # 軌道計算に使う要素のみをGlobal変数として抽出    ##open-end 内で定義した変数はopen-endの外では使えないため
    global Epoch = parse(Float64, line1[4])
    print("Epoch:元期 : ")
    println(Epoch)
    global Incl = parse(Float64, line2[3])
    print("Incliation:軌道傾斜角 : ")
    println(Epoch)
    global RAAN = parse(Float64, line2[4])
    print("Right Ascension of Ascending Node:昇交点赤経 : ")
    println(RAAN)
    global Ecce = parse(Float64, line2[5])/(10^length(line2[5]))
    print("Eccentricity:離心率 : ")
    println(Ecce)
    global AP = parse(Float64, line2[6])
    print("Argument Perigee:近地点離隔 : ")
    println(AP)
    global MA = parse(Float64, line2[7])
    print("Mean Anomaly:平均近点角 : ")
    println(MA)
    # 通算周回数が増え10000回を超えると2行目の8要素と9要素が結合するため、例外処理
    if(size(line2)[1]<9)
        global MM = parse(Float64, line2[8][1:length(line2[8])-6])
        print("Mean Motion:平均運動 : ")
        println(MM)
        global RNE = parse(Int32, line2[8][length(line2[8])-5:length(line2[8])-1])
        print("Revolution Number at Epoch:通算周回数 : ")
        println(RNE)
    else
        global MM = parse(Float64, line2[8])
        print("Mean Motion:平均運動 : ")
        println(MM)
        global RevNum = parse(Int32, line2[9][1:length(line2[9])-1])
        print("Revolution Number at Epoch:通算周回数 : ")
        println(RNE)
    end
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim