text:
###
### Define data structures used by later tests.
###

struct UnknownScalarSet{T} <: MOI.AbstractScalarSet
    constant::T
end

MOI.constant(set::UnknownScalarSet) = set.constant

Base.copy(set::UnknownScalarSet) = UnknownScalarSet(copy(MOI.constant(set)))

MOIU.supports_shift_constant(::Type{<:UnknownScalarSet}) = true
function MOIU.shift_constant(set::UnknownScalarSet, offset)
    return UnknownScalarSet(MOI.constant(set) + offset)
end

struct UnknownVectorSet <: MOI.AbstractVectorSet end

abstract type BadModel{T} <: MOI.ModelLike end

function MOI.get(::BadModel, ::MOI.ListOfModelAttributesSet)
    return MOI.AbstractModelAttribute[]
end

MOI.get(::BadModel, ::MOI.ListOfVariableIndices) = [MOI.VariableIndex(1)]

function MOI.get(::BadModel, ::MOI.ListOfVariableAttributesSet)
    return MOI.AbstractVariableAttribute[]
end

function MOI.get(::BadModel, ::MOI.ListOfConstraintAttributesSet)
    return MOI.AbstractConstraintAttribute[]
end

function MOI.get(::BadModel{T}, ::MOI.ListOfConstraintTypesPresent) where {T}
    return [(MOI.VariableIndex, MOI.EqualTo{T})]
end

function MOI.get(::BadModel, ::MOI.ListOfConstraintIndices{F,S}) where {F,S}
    return [MOI.ConstraintIndex{F,S}(1)]
end

function MOI.get(
    ::BadModel{T},
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}},
) where {T}
    return MOI.VariableIndex(1)
end

function MOI.get(
    ::BadModel{T},
    ::MOI.ConstraintSet,
    ::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}},
) where {T}
    return MOI.EqualTo(T(0))
end

struct BadConstraintModel{T} <: BadModel{T}
    BadConstraintModel(T = Float64) = new{T}()
end

function MOI.get(
    ::BadConstraintModel{T},
    ::MOI.ListOfConstraintTypesPresent,
) where {T}
    return [
        (MOI.VariableIndex, MOI.EqualTo{T}),
        (MOI.VariableIndex, UnknownScalarSet{T}),
    ]
end

function MOI.get(
    ::BadModel,
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.VariableIndex,<:UnknownScalarSet},
)
    return MOI.VariableIndex(1)
end

function MOI.get(
    ::BadModel,
    ::MOI.ConstraintSet,
    ::MOI.ConstraintIndex{MOI.VariableIndex,UnknownScalarSet{T}},
) where {T}
    return UnknownScalarSet(T(1))
end

struct UnknownModelAttribute <: MOI.AbstractModelAttribute end

struct BadModelAttributeModel{T} <: BadModel{T}
    BadModelAttributeModel(T = Float64) = new{T}()
end

MOI.get(src::BadModelAttributeModel, ::UnknownModelAttribute) = 0

function MOI.get(::BadModelAttributeModel, ::MOI.ListOfModelAttributesSet)
    return MOI.AbstractModelAttribute[UnknownModelAttribute()]
end

struct UnknownVariableAttribute <: MOI.AbstractVariableAttribute end
struct BadVariableAttributeModel{T} <: BadModel{T}
    BadVariableAttributeModel(T = Float64) = new{T}()
end

function MOI.get(
    ::BadVariableAttributeModel,
    ::UnknownVariableAttribute,
    ::MOI.VariableIndex,
)
    return 0
end

function MOI.get(::BadVariableAttributeModel, ::MOI.ListOfVariableAttributesSet)
    return MOI.AbstractVariableAttribute[UnknownVariableAttribute()]
end

struct UnknownConstraintAttribute <: MOI.AbstractConstraintAttribute end

struct BadConstraintAttributeModel{T} <: BadModel{T}
    BadConstraintAttributeModel(T = Float64) = new{T}()
end

function MOI.get(
    ::BadConstraintAttributeModel,
    ::UnknownConstraintAttribute,
    ::MOI.ConstraintIndex,
)
    return 0
end

function MOI.get(
    ::BadConstraintAttributeModel,
    ::MOI.ListOfConstraintAttributesSet,
)
    return MOI.AbstractConstraintAttribute[UnknownConstraintAttribute()]
end

###
### The actual tests
###

"""
    test_model_default_ObjectiveSense(model::MOI.ModelLike, ::Config)

Test that the default ObjectiveSense is FEASIBILITY_SENSE.
"""
function test_model_default_ObjectiveSense(model::MOI.ModelLike, ::Config)
    MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    _test_attribute_value_type(model, MOI.ObjectiveSense())
    return
end

"""
    test_model_default_TerminationStatus(model::MOI.AbstractOptimizer, ::Config)

Test that the default TerminationStatus is OPTIMIZE_NOT_CALLED.
"""
function test_model_default_TerminationStatus(
    model::MOI.AbstractOptimizer,
    ::Config,
)
    MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    _test_attribute_value_type(model, MOI.TerminationStatus())
    return
end

test_model_default_TerminationStatus(::MOI.ModelLike, ::Config) = nothing

"""
    test_model_default_PrimalStatus(model::MOI.AbstractOptimizer, ::Config)

Test that the default PrimalStatus is NO_SOLUTION.
"""
function test_model_default_PrimalStatus(model::MOI.AbstractOptimizer, ::Config)
    MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    _test_attribute_value_type(model, MOI.PrimalStatus())
    return
end

test_model_default_PrimalStatus(::MOI.ModelLike, ::Config) = nothing

"""
    test_model_default_DualStatus(model::MOI.AbstractOptimizer, ::Config)

Test that the default DualStatus is NO_SOLUTION.
"""
function test_model_default_DualStatus(model::MOI.AbstractOptimizer, ::Config)
    MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION
    _test_attribute_value_type(model, MOI.DualStatus())
    return
end

test_model_default_DualStatus(::MOI.ModelLike, ::Config) = nothing

"""
    test_model_VariableName(model::MOI.ModelLike, ::Config)

Test MOI.VariableName.
"""
function test_model_VariableName(model::MOI.ModelLike, ::Config)
    @requires MOI.supports(model, MOI.VariableName(), MOI.VariableIndex)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.VariableName(), x[1], "x1")
    @test MOI.get(model, MOI.VariableIndex, "x1") == x[1]
    MOI.set(model, MOI.VariableName(), x[1], "x2")
    @test MOI.get(model, MOI.VariableIndex, "x1") === nothing
    @test MOI.get(model, MOI.VariableIndex, "x2") == x[1]
    MOI.set(model, MOI.VariableName(), x[2], "x1")
    @test MOI.get(model, MOI.VariableIndex, "x1") == x[2]
    MOI.set(model, MOI.VariableName(), x[1], "x1")
    @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x1")
    _test_attribute_value_type(model, MOI.VariableName(), x[1])
    return
end

"""
    test_model_VariableIndex_ConstraintName(
        model::MOI.ModelLike,
        ::Config,
    )

Test ConstraintName for VariableIndex constraints.
"""
function test_model_VariableIndex_ConstraintName(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, x, MOI.GreaterThan(T(0)))
    @test_throws(
        MOI.VariableIndexConstraintNameError(),
        MOI.set(model, MOI.ConstraintName(), c, "c1"),
    )
    return
end

"""
    test_model_ScalarAffineFunction_ConstraintName(
        model::MOI.ModelLike,
        ::Config,
    )

Test ConstraintName for ScalarAffineFunction constraints.
"""
function test_model_ScalarAffineFunction_ConstraintName(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0))
    c1 = MOI.add_constraint(model, f, MOI.GreaterThan(T(0)))
    c2 = MOI.add_constraint(model, f, MOI.LessThan(T(1)))
    MOI.set(model, MOI.ConstraintName(), c1, "c1")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") == c1
    MOI.set(model, MOI.ConstraintName(), c1, "c2")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") === nothing
    @test MOI.get(model, MOI.ConstraintIndex, "c2") == c1
    MOI.set(model, MOI.ConstraintName(), c2, "c1")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") == c2
    MOI.set(model, MOI.ConstraintName(), c1, "c1")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "c1")
    return
end

"""
    test_model_Name(model::MOI.ModelLike, ::Config)

Test the `MOI.Name` attribute.
"""
function test_model_Name(model::MOI.ModelLike, ::Config)
    @requires MOI.supports(model, MOI.Name())
    @test !(MOI.Name() in MOI.get(model, MOI.ListOfModelAttributesSet()))
    @test MOI.get(model, MOI.Name()) == ""
    MOI.set(model, MOI.Name(), "Name1")
    @test MOI.Name() in MOI.get(model, MOI.ListOfModelAttributesSet())
    @test MOI.get(model, MOI.Name()) == "Name1"
    MOI.set(model, MOI.Name(), "Name2")
    @test MOI.Name() in MOI.get(model, MOI.ListOfModelAttributesSet())
    @test MOI.get(model, MOI.Name()) == "Name2"
    _test_attribute_value_type(model, MOI.Name())
    return
end

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

Test a variety of name attribute behavior.
"""
function test_model_Name_VariableName_ConstraintName(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_incremental_interface(model)
    @requires MOI.supports(model, MOI.VariableName(), MOI.VariableIndex)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    @test MOI.get(
        model,
        MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.LessThan{T}}(),
    ) == 0
    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.VariableName(), v[1]) == ""
    x, cx = MOI.add_constrained_variable(model, MOI.GreaterThan(T(0)))
    @test MOI.get(model, MOI.VariableName(), x) == ""
    y, cy = MOI.add_constrained_variables(model, MOI.Nonpositives(4))
    for yi in y
        @test MOI.get(model, MOI.VariableName(), yi) == ""
    end
    @test MOI.get(model, MOI.ConstraintName(), cy) == ""
    MOI.set(model, MOI.VariableName(), v[1], "")
    MOI.set(model, MOI.VariableName(), v[2], "") # Shouldn't error with duplicate empty name
    MOI.set(model, MOI.VariableName(), x, "")
    for yi in y
        MOI.set(model, MOI.VariableName(), yi, "")
    end
    MOI.set(model, MOI.VariableName(), v[1], "Var1")
    MOI.set(model, MOI.VariableName(), v[2], "Var1")
    # Lookup must fail when there are multiple variables with the same name.
    @test_throws Exception MOI.get(model, MOI.VariableIndex, "Var1")
    MOI.set(model, MOI.VariableName(), v[2], "Var2")
    @test MOI.get(model, MOI.VariableIndex, "Var1") == v[1]
    @test MOI.get(model, MOI.VariableIndex, "Var2") == v[2]
    @test MOI.get(model, MOI.VariableIndex, "Var3") === nothing
    MOI.set(model, MOI.VariableName(), x, "Var1")
    @test_throws Exception MOI.get(model, MOI.VariableIndex, "Var1")
    MOI.set(model, MOI.VariableName(), x, "Varx")
    @test MOI.get(model, MOI.VariableIndex, "Var1") == v[1]
    @test MOI.get(model, MOI.VariableIndex, "Var2") == v[2]
    @test MOI.get(model, MOI.VariableIndex, "Varx") == x
    @test MOI.get(model, MOI.VariableIndex, "Var3") === nothing
    vynames = ["VarX", "Var2", "Vary1", "Vary2", "Vary3", "Vary4"]
    MOI.set(model, MOI.VariableName(), [v; y], vynames)
    @test MOI.get(model, MOI.VariableName(), v) == vynames[1:2]
    @test MOI.get(model, MOI.VariableName(), y) == vynames[3:6]
    @test MOI.get(model, MOI.VariableName(), [v; y]) == vynames
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(T(1), v), T(0)),
        MOI.LessThan(T(1)),
    )
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    c2 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(T[-1, 1], v), T(0)),
        MOI.EqualTo(T(0)),
    )
    @test MOI.get(model, MOI.ConstraintName(), c) == ""
    @test MOI.get(model, MOI.ConstraintName(), c2) == ""
    @test MOI.get(model, MOI.ConstraintName(), cy) == ""
    @requires MOI.supports(model, MOI.ConstraintName(), typeof(c))
    MOI.set(model, MOI.ConstraintName(), c, "")
    @requires MOI.supports(model, MOI.ConstraintName(), typeof(c2))
    MOI.set(model, MOI.ConstraintName(), c2, "") # Shouldn't error with duplicate empty name
    @requires MOI.supports(model, MOI.ConstraintName(), typeof(cy))
    MOI.set(model, MOI.ConstraintName(), cy, "")
    MOI.set(model, MOI.ConstraintName(), c, "Con0")
    @test MOI.get(model, MOI.ConstraintName(), c) == "Con0"
    MOI.set(model, MOI.ConstraintName(), c2, "Con0")
    # Lookup must fail when multiple constraints have the same name.
    @test_throws Exception MOI.get(model, MOI.ConstraintIndex, "Con0")
    @test_throws Exception MOI.get(model, typeof(c), "Con0")
    MOI.set(model, MOI.ConstraintName(), [c], ["Con1"])
    @test MOI.get(model, MOI.ConstraintName(), [c]) == ["Con1"]
    @test MOI.get(
        model,
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.LessThan{T}},
        "Con1",
    ) == c
    @test MOI.get(
        model,
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
        "Con1",
    ) === nothing
    @test MOI.get(
        model,
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}},
        "Con1",
    ) === nothing
    @test MOI.get(
        model,
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.LessThan{T}},
        "Con2",
    ) === nothing
    @test MOI.get(model, MOI.ConstraintIndex, "Con1") == c
    @test MOI.get(model, MOI.ConstraintIndex, "Con2") === nothing
    MOI.set(model, MOI.ConstraintName(), [c2, cy], ["Con2", "Con2"])
    @test_throws Exception MOI.get(model, MOI.ConstraintIndex, "Con2")
    @test_throws Exception MOI.get(model, typeof(c2), "Con2")
    @test_throws Exception MOI.get(model, typeof(cy), "Con2")
    MOI.set(model, MOI.ConstraintName(), cy, "Con4")
    for (i, ca) in zip([1, 2, 4], [c, c2, cy])
        namea = "Con$i"
        @test MOI.get(model, MOI.ConstraintName(), ca) == namea
        @test MOI.get(model, typeof(ca), namea) == ca
        @test MOI.get(model, MOI.ConstraintIndex, namea) == ca
        for cb in [c, c2, cy]
            if ca === cb
                continue
            end
            nameb = MOI.get(model, MOI.ConstraintName(), cb)
            @test MOI.get(model, typeof(cb), namea) === nothing
            @test MOI.get(model, typeof(ca), nameb) === nothing
        end
    end
    if _supports(config, MOI.delete)
        MOI.delete(model, v[2])
        @test MOI.get(model, MOI.VariableIndex, "Var2") === nothing
        MOI.delete(model, c)
        @test MOI.get(model, typeof(c), "Con1") === nothing
        @test MOI.get(model, MOI.ConstraintIndex, "Con1") === nothing
        MOI.delete(model, x)
        @test MOI.get(model, MOI.VariableIndex, "Varx") === nothing
        @test MOI.get(model, MOI.ConstraintIndex, "Con3") === nothing
        @test MOI.get(model, typeof(c2), "Con2") === c2
        @test MOI.get(model, MOI.ConstraintIndex, "Con2") === c2
        MOI.delete(model, y)
        @test MOI.get(model, typeof(cy), "Con4") === nothing
        @test MOI.get(model, MOI.ConstraintIndex, "Con4") === nothing
        for i in 1:4
            @test MOI.get(model, MOI.VariableIndex, "Vary$i") === nothing
        end
        MOI.set(model, MOI.ConstraintName(), c2, "Con4")
        @test MOI.get(model, typeof(c2), "Con4") === c2
        @test MOI.get(model, MOI.ConstraintIndex, "Con4") === c2
    end
    return
end

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

Test duplicate variable names.
"""
function test_model_duplicate_VariableName(model::MOI.ModelLike, config::Config)
    (x, y, z) = MOI.add_variables(model, 3)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.VariableName(), y, "x")
    MOI.set(model, MOI.VariableName(), z, "z")
    @test MOI.get(model, MOI.VariableIndex, "z") == z
    @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x")
    MOI.set(model, MOI.VariableName(), y, "y")
    @test MOI.get(model, MOI.VariableIndex, "x") == x
    @test MOI.get(model, MOI.VariableIndex, "y") == y
    MOI.set(model, MOI.VariableName(), z, "x")
    @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x")
    if _supports(config, MOI.delete)
        MOI.delete(model, x)
        @test MOI.get(model, MOI.VariableIndex, "x") == z
    end
    return
end

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

Test duplicate names in ScalarAffineFunction constraints.
"""
function test_model_duplicate_ScalarAffineFunction_ConstraintName(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variables(model, 3)
    fs = [
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), xi)], T(0)) for
        xi in x
    ]
    c = MOI.add_constraints(model, fs, MOI.GreaterThan(T(0)))
    MOI.set(model, MOI.ConstraintName(), c[1], "x")
    MOI.set(model, MOI.ConstraintName(), c[2], "x")
    MOI.set(model, MOI.ConstraintName(), c[3], "z")
    @test MOI.get(model, MOI.ConstraintIndex, "z") == c[3]
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
    MOI.set(model, MOI.ConstraintName(), c[2], "y")
    @test MOI.get(model, MOI.ConstraintIndex, "x") == c[1]
    @test MOI.get(model, MOI.ConstraintIndex, "y") == c[2]
    MOI.set(model, MOI.ConstraintName(), c[3], "x")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
    if _supports(config, MOI.delete)
        MOI.delete(model, c[1])
        @test MOI.get(model, MOI.ConstraintIndex, "x") == c[3]
    end
    return
end

"""
    test_model_is_valid(model::MOI.ModelLike, config::Config{T}) where {T}

Test various parts of `MOI.is_valid`.

Taken from https://github.com/jump-dev/MathOptInterfaceUtilities.jl/issues/41
"""
function test_model_is_valid(model::MOI.ModelLike, config::Config{T}) where {T}
    @requires MOI.supports_incremental_interface(model)
    v = MOI.add_variables(model, 2)
    @test MOI.is_valid(model, v[1])
    @test MOI.is_valid(model, v[2])
    x = MOI.add_variable(model)
    @test MOI.is_valid(model, x)
    if _supports(config, MOI.delete)
        MOI.delete(model, x)
        @test !MOI.is_valid(model, x)
    end
    cf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(T(1), v), T(0))
    @requires MOI.supports_constraint(model, typeof(cf), MOI.LessThan{T})
    c = MOI.add_constraint(model, cf, MOI.LessThan(T(1)))
    @test MOI.is_valid(model, c)
    SAF = MOI.ScalarAffineFunction
    @test !MOI.is_valid(
        model,
        MOI.ConstraintIndex{SAF{Float32},MOI.LessThan{Float32}}(1),
    )
    @test !MOI.is_valid(
        model,
        MOI.ConstraintIndex{SAF{Float32},MOI.LessThan{T}}(1),
    )
    @test !MOI.is_valid(
        model,
        MOI.ConstraintIndex{SAF{T},MOI.LessThan{Float32}}(1),
    )
    VQF = MOI.VectorQuadraticFunction
    @test !MOI.is_valid(
        model,
        MOI.ConstraintIndex{VQF{T},MOI.SecondOrderCone}(1),
    )
    return
end

"""
    test_model_empty(model::MOI.ModelLike, ::Config)

Test `MOI.empty` and `MOI.is_empty!`.
"""
function test_model_empty(model::MOI.ModelLike, ::Config)
    @test MOI.is_empty(model)
    x = MOI.add_variable(model)
    @test !MOI.is_empty(model)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    return
end

"""
    test_model_copy_to_UnsupportedConstraint(model::MOI.ModelLike, ::Config)

Test an error is thrown when a constraint is unsupported.

!!! note
    This function uses the two-argument `MOI.optimize!`. This ensures that we
    test three cases:
     * models implementing the two-argument method
     * models using the generic fallback and erroring on `copy_to`
     * models using the generic fallback and erroring on `optimize!`.
"""
function test_model_copy_to_UnsupportedConstraint(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @test !MOI.supports_constraint(
        model,
        MOI.VariableIndex,
        UnknownScalarSet{T},
    )
    @test_throws(
        MOI.UnsupportedConstraint,
        MOI.optimize!(model, BadConstraintModel()),
    )
    return
end

"""
    test_model_copy_to_UnsupportedAttribute(model::MOI.ModelLike, ::Config)

Test an error is thrown when an attribute is unsupported.

!!! note
    This function uses the two-argument `MOI.optimize!`. This ensures that we
    test three cases:
     * models implementing the two-argument method
     * models using the generic fallback and erroring on `copy_to`
     * models using the generic fallback and erroring on `optimize!`.
"""
function test_model_copy_to_UnsupportedAttribute(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    # ModelAttribute
    @test !MOI.supports(model, UnknownModelAttribute())
    @test_throws(
        MOI.UnsupportedAttribute,
        MOI.optimize!(model, BadModelAttributeModel(T)),
    )
    # VariableAttribute
    @test !MOI.supports(model, UnknownVariableAttribute(), MOI.VariableIndex)
    @test_throws(
        MOI.UnsupportedAttribute,
        MOI.optimize!(model, BadVariableAttributeModel(T)),
    )
    # ConstraintAttribute
    @test !MOI.supports(
        model,
        UnknownConstraintAttribute(),
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}},
    )
    @test_throws(
        MOI.UnsupportedAttribute,
        MOI.optimize!(model, BadConstraintAttributeModel(T)),
    )
    return
end

function test_model_supports_constraint_VariableIndex_EqualTo(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{T})
    # Pick a "bad" coefficient type that should fail tests.
    @test !MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{UInt8})
    # Scalar-in-vector
    @test !MOI.supports_constraint(model, MOI.VariableIndex, MOI.Zeros)
    return
end

function test_model_supports_constraint_ScalarAffineFunction_EqualTo(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    # Pick a "bad" coefficient type that should fail tests.
    @test !MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{UInt8},
        MOI.EqualTo{UInt8},
    )
    return
end

function test_model_supports_constraint_VectorOfVariables_Nonnegatives(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.Nonnegatives,
    )
    # Pick a "bad" coefficient type that should fail tests.
    @test !MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        UnknownVectorSet,
    )
    @test !MOI.supports_constraint(model, MOI.VectorOfVariables, MOI.EqualTo{T})
    return
end

"""
    test_model_ordered_indices(model::MOI.ModelLike, ::Config{T}) where {T}

Test whether the model returns ListOfVariableIndices and ListOfConstraintIndices
sorted by creation time.
"""
function test_model_ordered_indices(model::MOI.ModelLike, ::Config{T}) where {T}
    @requires MOI.supports_incremental_interface(model)
    v1 = MOI.add_variable(model)
    @test MOI.get(model, MOI.ListOfVariableIndices()) == [v1]
    v2 = MOI.add_variable(model)
    @test MOI.get(model, MOI.ListOfVariableIndices()) == [v1, v2]
    MOI.delete(model, v1)
    @test MOI.get(model, MOI.ListOfVariableIndices()) == [v2]
    v3 = MOI.add_variable(model)
    @test MOI.get(model, MOI.ListOfVariableIndices()) == [v2, v3]
    v4 = MOI.add_variable(model)
    @test MOI.get(model, MOI.ListOfVariableIndices()) == [v2, v3, v4]
    # Note: there are too many combinations to test, so we're just going to
    # check VariableIndex-in-LessThan and hope it works for the rest
    c1 = MOI.add_constraint(model, v2, MOI.LessThan(T(1)))
    @test c1.value == v2.value
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}(),
    ) == [c1]
    c2 = MOI.add_constraint(model, v3, MOI.LessThan(T(2)))
    @test c2.value == v3.value
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}(),
    ) == [c1, c2]
    MOI.delete(model, c1)
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}(),
    ) == [c2]
    c3 = MOI.add_constraint(model, v4, MOI.LessThan(T(3)))
    @test c3.value == v4.value
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}(),
    ) == [c2, c3]
    return
end

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

Test adding a linear constraint with a non-zero function constant.

This should either work, or error with `MOI.ScalarFunctionConstantNotZero` if
the model does not support it.

To skip this test, pass `MOI.ScalarFunctionConstantNotZero` to the `exclude`
argument of [`Config`](@ref).
"""
function test_model_ScalarFunctionConstantNotZero(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires _supports(config, MOI.ScalarFunctionConstantNotZero)
    function _error(S, value)
        F = MOI.ScalarAffineFunction{T}
        return MOI.ScalarFunctionConstantNotZero{T,F,S}(value)
    end
    try
        f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], T(1))
        c = MOI.add_constraint(model, f, MOI.EqualTo(T(2)))
        @requires _supports(config, MOI.ConstraintFunction)
        @test MOI.get(model, MOI.ConstraintFunction(), c) ≈ f
    catch err
        @test err == _error(MOI.EqualTo{T}, T(1))
    end
    try
        f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], T(2))
        c = MOI.add_constraint(model, f, MOI.GreaterThan(T(1)))
        @requires _supports(config, MOI.ConstraintFunction)
        @test MOI.get(model, MOI.ConstraintFunction(), c) ≈ f
    catch err
        @test err == _error(MOI.GreaterThan{T}, T(2))
    end
    return
end

"""
    test_model_LowerBoundAlreadySet(
        model::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Test that setting a lower-bound twice throws `LowerBoundAlreadySet`.
"""
function test_model_LowerBoundAlreadySet(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.VariableIndex,
        MOI.GreaterThan{T},
    )
    @requires _supports(config, MOI.delete)
    x = MOI.add_variable(model)
    lb = T(0)
    sets = [MOI.EqualTo(lb), MOI.Interval(lb, lb)]
    set2 = MOI.GreaterThan(lb)
    for set1 in sets
        if !MOI.supports_constraint(model, MOI.VariableIndex, typeof(set1))
            continue
        end
        ci = MOI.add_constraint(model, x, set1)
        err = MOI.LowerBoundAlreadySet{typeof(set1),typeof(set2)}(x)
        @test_throws err MOI.add_constraint(model, x, set2)
        MOI.delete(model, ci)
        ci = MOI.add_constraint(model, x, set2)
        err = MOI.LowerBoundAlreadySet{typeof(set2),typeof(set1)}(x)
        @test_throws err MOI.add_constraint(model, x, set1)
        MOI.delete(model, ci)
    end
    return
end

"""
    test_model_UpperBoundAlreadySet(
        model::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Test that setting a lower-bound twice throws `UpperBoundAlreadySet`.
"""
function test_model_UpperBoundAlreadySet(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    x = MOI.add_variable(model)
    ub = T(0)
    @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.LessThan{T})
    sets = [MOI.EqualTo(ub), MOI.Interval(ub, ub)]
    set2 = MOI.LessThan(ub)
    for set1 in sets
        if !MOI.supports_constraint(model, MOI.VariableIndex, typeof(set1))
            continue
        end
        ci = MOI.add_constraint(model, x, set1)
        err = MOI.UpperBoundAlreadySet{typeof(set1),typeof(set2)}(x)
        @test_throws err MOI.add_constraint(model, x, set2)
        MOI.delete(model, ci)
        ci = MOI.add_constraint(model, x, set2)
        err = MOI.UpperBoundAlreadySet{typeof(set2),typeof(set1)}(x)
        @test_throws err MOI.add_constraint(model, x, set1)
        MOI.delete(model, ci)
    end
    return
end

"""
    test_model_delete(model::MOI.ModelLike, config::Config{T}) where {T}

Test various operations to do with deleting variables and constraints.
"""
function test_model_delete(model::MOI.ModelLike, config::Config{T}) where {T}
    @requires _supports(config, MOI.delete)
    x = MOI.add_variable(model)
    cx = MOI.add_constraint(model, x, MOI.GreaterThan(T(0)))
    y = MOI.add_variables(model, 4)
    cy = MOI.add_constraint(model, y, MOI.Nonpositives(4))
    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, y[1])
    @test MOI.is_valid(model, y[2])
    @test MOI.is_valid(model, y[3])
    @test MOI.is_valid(model, y[4])
    @test MOI.is_valid(model, cx)
    @test MOI.is_valid(model, cy)
    @test MOI.get(model, MOI.ConstraintFunction(), cx) == x
    @test MOI.get(model, MOI.ConstraintSet(), cx) == MOI.GreaterThan(T(0))
    @test MOI.get(model, MOI.ConstraintFunction(), cy) ==
          MOI.VectorOfVariables(y)
    @test MOI.get(model, MOI.ConstraintSet(), cy) == MOI.Nonpositives(4)
    @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
        (MOI.VariableIndex, MOI.GreaterThan{T}),
        (MOI.VectorOfVariables, MOI.Nonpositives),
    ])
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(),
    ) == [cx]
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonpositives}(),
    ) == [cy]
    MOI.delete(model, y[3])
    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, y[1])
    @test MOI.is_valid(model, y[2])
    @test !MOI.is_valid(model, y[3])
    @test MOI.is_valid(model, y[4])
    @test MOI.is_valid(model, cx)
    @test MOI.is_valid(model, cy)
    @test MOI.get(model, MOI.ConstraintFunction(), cx) == x
    @test MOI.get(model, MOI.ConstraintSet(), cx) == MOI.GreaterThan(T(0))
    @test MOI.get(model, MOI.ConstraintFunction(), cy) ==
          MOI.VectorOfVariables(y[[1, 2, 4]])
    @test MOI.get(model, MOI.ConstraintSet(), cy) == MOI.Nonpositives(3)
    @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
        (MOI.VariableIndex, MOI.GreaterThan{T}),
        (MOI.VectorOfVariables, MOI.Nonpositives),
    ])
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(),
    ) == [cx]
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonpositives}(),
    ) == [cy]
    MOI.delete(model, y[1])
    @test MOI.is_valid(model, x)
    @test !MOI.is_valid(model, y[1])
    @test MOI.is_valid(model, y[2])
    @test !MOI.is_valid(model, y[3])
    @test MOI.is_valid(model, y[4])
    @test MOI.is_valid(model, cx)
    @test MOI.is_valid(model, cy)
    @test MOI.get(model, MOI.ConstraintFunction(), cx) == x
    @test MOI.get(model, MOI.ConstraintSet(), cx) == MOI.GreaterThan(T(0))
    @test MOI.get(model, MOI.ConstraintFunction(), cy) ==
          MOI.VectorOfVariables(y[[2, 4]])
    @test MOI.get(model, MOI.ConstraintSet(), cy) == MOI.Nonpositives(2)
    @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
        (MOI.VariableIndex, MOI.GreaterThan{T}),
        (MOI.VectorOfVariables, MOI.Nonpositives),
    ])
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(),
    ) == [cx]
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonpositives}(),
    ) == [cy]
    MOI.delete(model, x)
    @test !MOI.is_valid(model, x)
    @test !MOI.is_valid(model, y[1])
    @test MOI.is_valid(model, y[2])
    @test !MOI.is_valid(model, y[3])
    @test MOI.is_valid(model, y[4])
    @test !MOI.is_valid(model, cx)
    @test MOI.is_valid(model, cy)
    @test MOI.get(model, MOI.ConstraintFunction(), cy) ==
          MOI.VectorOfVariables(y[[2, 4]])
    @test MOI.get(model, MOI.ConstraintSet(), cy) == MOI.Nonpositives(2)
    @test MOI.get(model, MOI.ListOfConstraintTypesPresent()) ==
          [(MOI.VectorOfVariables, MOI.Nonpositives)]
    @test isempty(
        MOI.get(
            model,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(),
        ),
    )
    @test MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonpositives}(),
    ) == [cy]
    MOI.delete(model, y[[2, 4]])
    @test !MOI.is_valid(model, x)
    @test !MOI.is_valid(model, y[1])
    @test !MOI.is_valid(model, y[2])
    @test !MOI.is_valid(model, y[3])
    @test !MOI.is_valid(model, y[4])
    @test !MOI.is_valid(model, cx)
    @test !MOI.is_valid(model, cy)
    @test isempty(MOI.get(model, MOI.ListOfConstraintTypesPresent()))
    @test isempty(
        MOI.get(
            model,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(),
        ),
    )
    @test isempty(
        MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonpositives,
            }(),
        ),
    )
end

"""
    test_model_ListOfConstraintAttributesSet(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test issue #1429: ConstraintName should only be included in
`ListOfConstraintAttributesSet` if it actually is set.
"""
function test_model_ListOfConstraintAttributesSet(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports(
        model,
        MOI.ConstraintName(),
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
    )
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(T(0)))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0)),
        MOI.EqualTo(T(1)),
    )
    MOI.set(model, MOI.ConstraintName(), c, "c")
    @test MOI.get(
        model,
        MOI.ListOfConstraintAttributesSet{
            MOI.ScalarAffineFunction{T},
            MOI.EqualTo{T},
        }(),
    ) == [MOI.ConstraintName()]
    @test MOI.get(
        model,
        MOI.ListOfConstraintAttributesSet{
            MOI.VariableIndex,
            MOI.GreaterThan{T},
        }(),
    ) == []
    return
end

"""
    test_model_ModelFilter_AbstractModelAttribute(
        src::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Tests `Utilties.ModelFilter` of `AbstractModelAttribute`.
"""
function test_model_ModelFilter_AbstractModelAttribute(
    src::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports(src, MOI.Name())
    MOI.set(src, MOI.Name(), "src")
    dest = MOI.Utilities.Model{T}()
    MOI.copy_to(dest, MOI.Utilities.ModelFilter(src) do item
        return item != MOI.Name()
    end)
    @test MOI.get(dest, MOI.Name()) == ""
    return
end

"""
    test_model_ModelFilter_AbstractVariableAttribute(
        src::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Tests `Utilties.ModelFilter` of `AbstractVariableAttribute`.
"""
function test_model_ModelFilter_AbstractVariableAttribute(
    src::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports(src, MOI.VariableName(), MOI.VariableIndex)
    @requires MOI.supports(src, MOI.VariablePrimalStart(), MOI.VariableIndex)
    x = MOI.add_variable(src)
    MOI.set(src, MOI.VariableName(), x, "x")
    MOI.set(src, MOI.VariablePrimalStart(), x, T(1))
    dest = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}())
    index_map = MOI.copy_to(dest, MOI.Utilities.ModelFilter(src) do item
        return item != MOI.VariableName()
    end)
    @test MOI.get(dest, MOI.VariableName(), index_map[x]) == ""
    @test MOI.get(dest, MOI.VariablePrimalStart(), index_map[x]) == T(1)
    return
end

"""
    test_model_ModelFilter_AbstractConstraintAttribute(
        src::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Tests `Utilties.ModelFilter` of `AbstractConstraintAttribute`.
"""
function test_model_ModelFilter_AbstractConstraintAttribute(
    src::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires(
        MOI.supports(
            src,
            MOI.ConstraintName(),
            MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives},
        ),
    )
    @requires(
        MOI.supports(
            src,
            MOI.ConstraintDualStart(),
            MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives},
        ),
    )
    @requires(
        MOI.supports_constraint(src, MOI.VectorOfVariables, MOI.Nonnegatives),
    )
    x = MOI.add_variable(src)
    c = MOI.add_constraint(src, MOI.VectorOfVariables([x]), MOI.Nonnegatives(1))
    MOI.set(src, MOI.ConstraintName(), c, "c")
    MOI.set(src, MOI.ConstraintDualStart(), c, T[1])
    dest = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}())
    index_map = MOI.copy_to(dest, MOI.Utilities.ModelFilter(src) do item
        if item isa MOI.AbstractConstraintAttribute
            return item != MOI.ConstraintName()
        end
        return true
    end)
    @test MOI.get(dest, MOI.ConstraintName(), index_map[c]) == ""
    @test MOI.get(dest, MOI.ConstraintDualStart(), index_map[c]) == T[1]
    return
end

"""
    test_model_ModelFilter_ListOfConstraintIndices(
        src::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Tests `Utilties.ModelFilter` of `ListOfConstraintIndices`.
"""
function test_model_ModelFilter_ListOfConstraintIndices(
    src::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires(
        MOI.supports_constraint(src, MOI.VariableIndex, MOI.GreaterThan{T}),
    )
    x = MOI.add_variables(src, 4)
    MOI.add_constraint.(src, x, MOI.GreaterThan(T(0)))
    dest = MOI.Utilities.Model{T}()
    index_map = MOI.copy_to(dest, MOI.Utilities.ModelFilter(src) do item
        if item isa MOI.ConstraintIndex
            return isodd(item.value)
        end
        return true
    end)
    @test MOI.get(dest, MOI.NumberOfVariables()) == 4
    @test MOI.get(
        dest,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.GreaterThan{T}}(),
    ) == 2
    for xi in x
        @test haskey(index_map, xi)
        ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{T}}(xi.value)
        @test haskey(index_map, ci) == isodd(xi.value)
    end
    return
end

"""
    test_model_ModelFilter_ListOfConstraintTypesPresent(
        src::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Tests `Utilties.ModelFilter` of `ListOfConstraintTypesPresent`.
"""
function test_model_ModelFilter_ListOfConstraintTypesPresent(
    src::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires(
        MOI.supports_constraint(src, MOI.VariableIndex, MOI.GreaterThan{T}),
    )
    @requires(MOI.supports_constraint(src, MOI.VariableIndex, MOI.LessThan{T}))
    x = MOI.add_variables(src, 4)
    MOI.add_constraint.(src, x, MOI.GreaterThan(T(0)))
    MOI.add_constraint.(src, x, MOI.LessThan(T(1)))
    dest = MOI.Utilities.Model{T}()
    _ = MOI.copy_to(
        dest,
        MOI.Utilities.ModelFilter(src) do item
            if item isa Tuple{Type,Type}
                return item == (MOI.VariableIndex, MOI.LessThan{T})
            end
            return true
        end,
    )
    @test MOI.get(dest, MOI.NumberOfVariables()) == 4
    attr = MOI.NumberOfConstraints{MOI.VariableIndex,MOI.GreaterThan{T}}()
    @test MOI.get(dest, attr) == 0
    attr = MOI.NumberOfConstraints{MOI.VariableIndex,MOI.LessThan{T}}()
    @test MOI.get(dest, attr) == 4
    return
end


text:
## (c) 2016--2018 David A. van Leeuwen
## tests for ../src/index.jl

print("getindex, ")
## getindex
## Test Integer indices up to 5 dimensions, as well as CartesianIndexes
for i in 1:7
	dims = fill(3, i)
	global n1 = @inferred NamedArray(rand(dims...))
	for i in CartesianIndices(tuple(dims...))
		@test n1[i.I...] == n1.array[i.I...]
		@test n1[i] == n1.array[i]
		@test n1[[string(ii) for ii in i.I]...] == n1.array[i]
		n1[i.I...] +=  1 ## test setindex
	end
end

include("init-namedarrays.jl")

@test n[:] == view(n, :) == n.array[:] ## values
println(n)
println(names(n[:],1))
@test names(n[:], 1)[2] == ("two", "a")
@test names(n[:], 1)[5] == ("one", "c")
@test dimnames(n[:], 1) == (:A, :B)

n1 = @inferred NamedArray(rand(10))
@test n1[:] == n1
@test [x for x in n] == [x for x in n.array]

## more indexing
@test convert(Array, n["one", :]) == n.array[1,:]
@test n[Not("two"), :].array == n.array[1:1,:]
@test names(n[Not("two"), :]) == names(n[1:1, :])
@test n[:, ["b", "d"]] == view(n, :, ["b", "d"]) == n[:, [2, 4]]

@test names(n[Not("one"), :], 1) == ["two"]
@test names(n[1:1, Not("a")], 2) == ["b", "c", "d"]
## indexing by pair
@test n[:B=>"a", :A=>"two"] == n[2, 1]
@test n[:A=>"one", :B=>"d"] == n[1, 4]
@test n[:B => "b"] == n[:, "b"]
@test n[:A => "two"] == n["two", :]

## https://github.com/nalimilan/FreqTables.jl/issues/10
bi = [false, true, false, true] ## Array{Bool}
for i in 1:2
    @test n[:, bi] == view(n, :, bi) == n.array[:, bi]
    @test names(n[:, bi], 2) == ["b", "d"]
    global bi = BitArray(bi)
end

m = @inferred copy(n)
print("setindex, ")
## setindex
m[1,1] = 0
m[2,:] = 1:4
m[:,"c"] .= -1
m[1,[2,3]] = [10,20]
m["one", 4] = 5

@test m.array == [0. 10 20 5; 1 2 -1 4]
m2 = copy(m)
for i in eachindex(m)
	m[i] = 2*m[i]
end
@test 2*(m2.array) == m.array

m[:B=>"c", :A=>"one"] = π
@test m[1,3] == Float64(π)
m[:] = n
@test m == n

m[:] = 1:8
@test m[:] == collect(1:8)

m = @inferred NamedArray(rand(Int, 10))
m[2:5] .= -1
m[6:8] = 2:4
m[[1,9,10]] = 0:2
@test m == [0, -1, -1, -1, -1, 2, 3, 4, 1, 2]

n[1,1] = π
@test n.array[1,1] == Float64(π)

ext: jl


text:
module ValidationHelper

using Genie
using SearchLight.Validation

export output_errors

function output_errors(m, field::Symbol) :: String
  Validation.has_errors_for(m, field) ? """$(Validation.errors_to_string(m, field, "
\n", upper_case_first = true))""" : ""
end

end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

import MLJFlux
import MLJ
import DataFrames
import Statistics
import Flux
using Random
using PyPlot

Random.seed!(11)

features, targets = MLJ.@load_boston
features = DataFrames.DataFrame(features)
@show size(features)
@show targets[1:3]
first(features, 3) |> MLJ.pretty

train, test = MLJ.partition(MLJ.eachindex(targets), 0.70, rng=52)

mutable struct MyNetworkBuilder <: MLJFlux.Builder
    n1::Int #Number of cells in the first hidden layer
    n2::Int #Number of cells in the second hidden layer
end

function MLJFlux.build(model::MyNetworkBuilder, input_dims, output_dims)
    layer1 = Flux.Dense(input_dims, model.n1)
    layer2 = Flux.Dense(model.n1, model.n2)
    layer3 = Flux.Dense(model.n2, output_dims)
    return Flux.Chain(layer1, layer2, layer3)
end

myregressor = MyNetworkBuilder(20, 10)

nnregressor = MLJFlux.NeuralNetworkRegressor(builder=myregressor, epochs=10)

mach = MLJ.machine(nnregressor, features, targets)

MLJ.fit!(mach, rows=train, verbosity=3)

preds = MLJ.predict(mach, features[test, :])

print(preds[1:5])

nnregressor.epochs = 15

MLJ.fit!(mach, rows=train, verbosity=3)

nnregressor.batch_size = 2
MLJ.fit!(mach, rows=train, verbosity=3)

# Inspecting out-of-sample loss as a function of epochs

r = MLJ.range(nnregressor, :epochs, lower=1, upper=30, scale=:log10)
curve = MLJ.learning_curve(nnregressor, features, targets,
                       range=r,
                       resampling=MLJ.Holdout(fraction_train=0.7),
                       measure=MLJ.l2)

figure(figsize=(8,6))

plt.plot(curve.parameter_values,
    curve.measurements)

yscale("log")
xlabel(curve.parameter_name)
ylabel("l2")

bs = MLJ.range(nnregressor, :batch_size, lower=1, upper=5)

tm = MLJ.TunedModel(model=nnregressor, ranges=[bs, ], measure=MLJ.l2)

m = MLJ.machine(tm, features, targets)

MLJ.fit!(m)

MLJ.fitted_params(m).best_model.batch_size

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -1.447605847981205 <= q <= 6.893417062244793)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.9197484844984549 + 0.6950857860734887*$(x[1]) + 0.2312849177869909*$(x[2]) + 0.027956483294022227*$(x[3]) + -0.9150118369982039*$(x[4])))) + 1/(1 + exp(-(0.925394406424513 + -0.0010850578907235509*$(x[1]) + -0.7627274772102703*$(x[2]) + -0.08863835978903944*$(x[3]) + 0.6322970885722761*$(x[4])))) + 1/(1 + exp(-(0.8777627162088257 + 0.18677359533052185*$(x[1]) + -0.4414865478782102*$(x[2]) + 0.1273485474015157*$(x[3]) + -0.06081575688773588*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    ModularProduct, modularproduct,
    CartesianProduct, cartesianproduct

# Modular product

struct ModularProductNode <: AbstractNode
    g::Int
    h::Int
end

struct ModularProductEdge <: UndirectedEdge
    hasedge::Bool
end

struct ModularProduct <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{ModularProductNode}
    edgeattrs::Vector{ModularProductEdge}
    cache::Dict{Symbol,Any}
end

"""
    modularproduct(G::OrderedGraph, H::OrderedGraph) -> ModularProduct

Return the modular product of graphs G and H.
"""
function modularproduct(G::OrderedGraph, H::OrderedGraph;
            nodematcher=(g,h)->true,
            edgefilter=(g1,g2,h1,h2)->hasedge(G,g1,g2)==hasedge(H,h1,h2),
            kwargs...)
    product = ModularProduct([], [], [], [], Dict())
    ndict = Dict{Int,Dict{Int,Int}}() # Ref to node indices of the product
    # Modular product nodes
    nattrs = ModularProductNode[]
    for g in 1:nodecount(G)
        ndict[g] = Dict{Int,Int}()
        for h in 1:nodecount(H)
            ndict[g][h] = addnode!(product, ModularProductNode(g, h))
        end
    end
    (nodecount(G) < 2 || nodecount(H) < 2) && return product
    # Modular product edges
    for (g1, g2) in combinations(nodecount(G))
        for (h1, h2) in combinations(nodecount(H))
            edgefilter(g1, g2, h1, h2) || continue
            if nodematcher(g1, h1) && nodematcher(g2, h2)
                addedge!(
                    product, ndict[g1][h1], ndict[g2][h2],
                    ModularProductEdge(hasedge(G, g1, g2))
                )
            end
            if nodematcher(g1, h2) && nodematcher(g2, h1)
                addedge!(
                    product, ndict[g1][h2], ndict[g2][h1],
                    ModularProductEdge(hasedge(G, g1, g2))
                )
            end
        end
    end
    return product
end

# Cartesian product

struct CartesianProductNode <: AbstractNode
    g::Int
    h::Int
end

struct CartesianProduct <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{CartesianProductNode}
    cache::Dict{Symbol,Any}
end

"""
    cartesianproduct(G::OrderedGraph, H::OrderedGraph) -> CartesianProduct

Return the cartesian product of graphs G and H.
"""
function cartesianproduct(G::OrderedGraph, H::OrderedGraph)
    product = CartesianProduct([], [], [], Dict())
    ndict = Dict{Int,Dict{Int,Int}}() # Ref to node indices of the product
    #Cartesian product nodes
    nattrs = CartesianProductNode[]
    for g in 1:nodecount(G)
        ndict[g] = Dict{Int,Int}()
        for h in 1:nodecount(H)
            ndict[g][h] = addnode!(product, CartesianProductNode(g, h))
        end
    end
    (nodecount(G) < 2 || nodecount(H) < 2) && return product
    # Cartesian product edges
    for i in 1:nodecount(G)
        for (u, v) in edgesiter(H)
            addedge!(product, ndict[i][u], ndict[i][v])
        end
    end
    for i in 1:nodecount(H)
        for (u, v) in edgesiter(G)
            addedge!(product, ndict[u][i], ndict[v][i])
        end
    end
    return product
end


text:
import DataFormatter as DF
using Test
using DataFrames

@testset "2 times series TimeseriesDataset" begin
    tsdf = DataFrame(:t=>rand(6),:y=>rand(6),:v=>[1,1,1,2,2,2])
    tsds = DF.TimeseriesDataset(tsdf,:v)
    ts_groups = DF.split_timeseries(tsds)
    @test length(ts_groups) == 2
end

@testset "1 times series TimeseriesDataset" begin
    tsdf = DataFrame(:t=>rand(6),:y=>rand(6),:v=>[1,1,1,1,1,1])
    tsds = DF.TimeseriesDataset(tsdf,:v)
    ts_groups = DF.split_timeseries(tsds)
    @test length(ts_groups) == 1
end

@testset "2 times series TimeseriesDataset with keys" begin
    tsdf = DataFrame(:t=>rand(6),:y=>rand(6),:v=>[1,1,1,1,1,1])
    keys = ["t","y","v"]
    keys2 = [:t,:y,:v]
    tsds = DF.TimeseriesDataset(tsdf, :v ; keys=keys)
    tsds = DF.TimeseriesDataset(tsdf, :v ; keys=keys2)

    ts_groups = DF.split_timeseries(tsds)
    @test length(ts_groups) == 1
end

@testset "colname not in keys" begin
    tsdf = DataFrame(:t=>rand(6),:y=>rand(6),:v=>[1,1,1,1,1,1])
    keys = ["t","y"]
    @test_throws ArgumentError DF.TimeseriesDataset(tsdf, :v ; keys=keys)
end

max_stars_repo_name: Leooop/DataFormatter

max_issues_repo_name: Leooop/DataFormatter

max_forks_repo_name: Leooop/DataFormatter


text:
doc"""
    split(x::Var, size::Vector)
    split(x::Var, dim::Int, size::Vector{Int})

# Example
```julia
T = Float32
x = Var(rand(T,10,10))
ys1 = split(x, [(5,10),(5,10)])
ys2 = split(x, 2, [2,3,5])
```
"""
function Base.split(x::Array{T,N}, size::Tuple) where {T,N}
    offset = 0
    map(size) do s
        p = pointer(x, offset+1)
        y = unsafe_wrap(typeof(x), p, s)
        y = Var(y, (split,x,offset))
        offset += length(y)
        y
    end
end

function Base.split(x::Var, dim::Int, size::Vector{Int})
    @assert sum(size) == Base.size(x,dim)
    if dim == ndims(x)
        offset = 0
        front = Base.front(Base.size(x))
        map(size) do s
            p = pointer(x.data, offset+1)
            y = unsafe_wrap(typeof(x.data), p, (front...,s))
            y = Var(y, (split,x,dim,offset))
            offset += length(y)
            y
        end
    else
        throw("Not implemented yet.")
    end
end
Base.split(x::Node, args...) = Node(split, x, args...)

function addgrad!(y::Var, ::typeof(split), x::Var, dim::Int, offset::Int)
    isvoid(x.grad) && return
    addto!(x.grad, offset, y.grad, 1, length(y))
end


text:
# ---
# title: 609. Find Duplicate File in System
# id: problem609
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, String
# link: <https://leetcode.com/problems/find-duplicate-file-in-system/description/>
# hidden: true
# ---
#
# Given a list of directory info including directory path, and all the files
# with contents in this directory, you need to find out all the groups of
# duplicate files in the file system in terms of their paths.
#
# A group of duplicate files consists of at least **two** files that have
# exactly the same content.
#
# A single directory info string in the **input** list has the following format:
#
# `"root/d1/d2/.../dm f1.txt(f1_content) f2.txt(f2_content) ...
# fn.txt(fn_content)"`
#
# It means there are **n** files (`f1.txt`, `f2.txt` ... `fn.txt` with content
# `f1_content`, `f2_content` ... `fn_content`, respectively) in directory
# `root/d1/d2/.../dm`. Note that n >= 1 and m >= 0. If m = 0, it means the
# directory is just the root directory.
#
# The **output** is a list of group of duplicate file paths. For each group, it
# contains all the file paths of the files that have the same content. A file
# path is a string that has the following format:
#
# `"directory_path/file_name.txt"`
#
# **Example 1:**
#
#
#
#     Input:
#     ["root/a 1.txt(abcd) 2.txt(efgh)", "root/c 3.txt(abcd)", "root/c/d 4.txt(efgh)", "root 4.txt(efgh)"]
#     Output:
#     [["root/a/2.txt","root/c/d/4.txt","root/4.txt"],["root/a/1.txt","root/c/3.txt"]]
#
#
#
#
# **Note:**
#
#   1. No order is required for the final output.
#   2. You may assume the directory name, file name and file content only has letters and digits, and the length of file content is in the range of [1,50].
#   3. The number of files given is in the range of [1,20000].
#   4. You may assume no files or directories share the same name in the same directory.
#   5. You may assume each given directory info represents a unique directory. Directory path and file info are separated by a single blank space.
#
#
#
# **Follow-up beyond contest:**
#
#   1. Imagine you are given a real file system, how will you search files? DFS or BFS?
#   2. If the file content is very large (GB level), how will you modify your solution?
#   3. If you can only read the file by 1kb each time, how will you modify your solution?
#   4. What is the time complexity of your modified solution? What is the most time-consuming part and memory consuming part of it? How to optimize?
#   5. How to make sure the duplicated files you find are not false positive?
#
#
## @lc code=start
using LeetCode

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


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using SpecialFunctions, Plots

logfact(n) = loggamma(n + 1)
stirling(n, ::Val{1}) = n*log(n)
stirling(n, ::Val{2}) = n*log(n) - n
stirling(n, ::Val{3}) = n*log(n) - n + (1/2)*log(n)
stirling(n, ::Val{4}) = n*log(n) - n + (1/2)*log(n) + log(√(2π))

PP = []
for (m, legend) in ((0, :bottomright), (1, :topleft), (2, :topleft), (23, :topleft))
    xmax = 10.0^m
    P = plot(; legend, xlim=(0, 1.05xmax), tickfontsize=6)
    plot!(logfact, 0.01, 1.05xmax; label="log n!")
    for k in 1:4
        plot!(x -> stirling(x, Val(k)), 0.01, 1.05xmax; label="stirling $k", ls=:auto)
    end
    push!(PP, P)
end
plot(PP...; size=(800, 600))

# %% [markdown]
# ```julia
# logfact(n) = loggamma(n + 1)
# stirling(n, ::Val{1}) = n*log(n)
# stirling(n, ::Val{2}) = n*log(n) - n
# stirling(n, ::Val{3}) = n*log(n) - n + (1/2)*log(n)
# stirling(n, ::Val{4}) = n*log(n) - n + (1/2)*log(n) + log(√(2π))
# ```

# %%

max_stars_repo_path: 0013/Stirling formula.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/Stirling formula.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/Stirling formula.jl

max_forks_repo_name: genkuroki/public


text:
module QXTns

using ITensors

include("utils.jl")
include("circuits.jl")
include("mock_tensor.jl")
include("tensor.jl")
include("tensor_network.jl")
include("tensor_network_circuit.jl")

end


text:
"""
This file have the basic outer approximation algoirthm
"""

"""
function solve_ref_nlp_model(model::MOI.AbstractOptimizer, oad::OAdata)
    This function solves reformulated nlp problem.
"""
function solve_ref_nlp_model(model::MOI.AbstractOptimizer, oad::OAdata)
    if oad.milp_sol_available
        vars = JuMP.all_variables(oad.ref_nlp_model)
        for i in 1:oad.ref_num_var
            if oad.ref_var_type[i] == :Int || oad.ref_var_type[i] ==:Bin
                JuMP.delete_lower_bound(vars[i])
                JuMP.delete_upper_bound(vars[i])
                JuMP.set_lower_bound(vars[i], oad.ref_mip_solution[i])
                JuMP.set_upper_bound(vars[i], oad.ref_mip_solution[i])
            end
        end
    end

    JuMP.set_optimizer(oad.ref_nlp_model, model.options.nlp_solver)
    JuMP.optimize!(oad.ref_nlp_model);
    #the problem is here
    if (JuMP.termination_status(oad.ref_nlp_model) == MOI.OPTIMAL) || (JuMP.termination_status(oad.ref_nlp_model) == MOI.LOCALLY_SOLVED)
        @info "\n The reformulated nlp problem has a feasible solution. \n"
        oad.milp_sol_available && (oad.ref_nlp_solution  = JuMP.value.(vars))
        !oad.milp_sol_available && (oad.ref_nlp_solution  = JuMP.value.(oad.ref_nlp_x))
        nlp_objval            = JuMP.objective_value(oad.ref_nlp_model)
        #update objval
        if (nlp_objval < oad.obj_val) && oad.milp_sol_available            #considering the problem is a min sense problem
            oad.obj_val       = nlp_objval
            oad.incumbent     = oad.ref_nlp_solution
            oad.new_incumbent = true
        end
        #update objgap
        if isfinite(oad.obj_val) && isfinite(oad.obj_bound)
            oad.obj_gap = (oad.obj_val - oad.obj_bound)/(abs(oad.obj_val) + 1e-5)
        end
    elseif (JuMP.termination_status(oad.ref_nlp_model) == MOI.INFEASIBLE) || (JuMP.termination_status(oad.ref_nlp_model) == MOI.DUAL_INFEASIBLE) || (JuMP.termination_status(oad.ref_nlp_model) == MOI.LOCALLY_INFEASIBLE) || (JuMP.termination_status(oad.ref_nlp_model) == MOI.INFEASIBLE_OR_UNBOUNDED)
        @info "\n The reformulated nlp problem is infeasible. \n"
        oad.nlp_infeasible = true
    end
end

"""
function solve_ref_mip_model(model::MOI.AbstractOptimizer, oad::OAdata)
    This function solves reformulatd mip problem.
"""
function solve_ref_mip_model(model::MOI.AbstractOptimizer, oad::OAdata)
    JuMP.set_optimizer(oad.ref_mip_model, model.options.mip_solver)
    JuMP.optimize!(oad.ref_mip_model);

    if (JuMP.termination_status(oad.ref_mip_model) ==  MOI.OPTIMAL) || (JuMP.termination_status(oad.ref_mip_model) ==  MOI.LOCALLY_SOLVED)
        @info "\n The reformulated mip problem has a feasible solution. \n"
        (oad.milp_sol_available == true) && (oad.prev_ref_mip_solution = oad.ref_mip_solution)
        (oad.milp_sol_available == false) && (oad.milp_sol_available = true)
        oad.ref_mip_solution = JuMP.value.(oad.ref_mip_x)
        # update bestbound
        mipobjbound = JuMP.objective_bound(oad.ref_mip_model)
        if isfinite(mipobjbound) && (mipobjbound > oad.obj_bound)
            oad.obj_bound = mipobjbound
        end
    else (JuMP.termination_status(oad.ref_mip_model) == MOI.INFEASIBLE) || (JuMP.termination_status(oad.ref_mip_model) == MOI.DUAL_INFEASIBLE) || (JuMP.termination_status(oad.ref_mip_model) == MOI.LOCALLY_INFEASIBLE) || (JuMP.termination_status(oad.ref_mip_model) == MOI.INFEASIBLE_OR_UNBOUNDED)
        @info "\n The reformulated mip problem is infeasible. \n"
        oad.mip_infeasible = true
    end
end

"""
function solve_ref_feasibility_model(model::MOI.AbstractOptimizer, oa_data::OAdata)
    This function solves reformulated feasibility problem.
"""
function solve_ref_feasibility_model(model::MOI.AbstractOptimizer, oad::OAdata)
    #fix int variables
    if oad.milp_sol_available
        for i in 1:oad.ref_num_var
            if (oad.ref_var_type[i] == :Int) || (oad.ref_var_type[i] ==:Bin)
                JuMP.set_lower_bound(oad.ref_feasibility_x[i], oad.ref_mip_solution[i])
                JuMP.set_upper_bound(oad.ref_feasibility_x[i], oad.ref_mip_solution[i])
            end
        end
    end
    #set optimizer and solve
    JuMP.set_optimizer(oad.ref_feasibility_model, model.options.nlp_solver)
    JuMP.optimize!(oad.ref_feasibility_model);

    if (JuMP.termination_status(oad.ref_feasibility_model) == MOI.OPTIMAL) || (JuMP.termination_status(oad.ref_feasibility_model) == MOI.LOCALLY_SOLVED)
        @info "\n The reformulated feasibility problem has a feasible solution. \n"
        oad.ref_feasibility_solution = JuMP.value.(oad.ref_feasibility_x)
    else
        @error "\n The reformulated feasibility problem is infeasible. \n"
    end
end

"""
function check_termination_condition(oad::OAdata)
    this function checks termination status of the oa algorithm.
"""
function check_termination_condition(oad::OAdata)
    # finish if optimal or cycling integer solutions
    if oad.obj_gap <= 0.00001
        oad.oa_status = :Optimal
    elseif round.(oad.prev_ref_mip_solution[oad.int_idx]) == round.(oad.ref_mip_solution[oad.int_idx])
        @warn "mixed-integer cycling detected, terminating outer approximation algorithm"
        if isfinite(oad.obj_gap)
            oad.oa_status = :Suboptimal
        else
            oad.oa_status = :FailedOA
        end
    elseif oad.mip_infeasible && isfinite(oad.obj_val)
        oad.oa_status = :Optimal
    end
end

"""
function oa_algorithm(model::MOI.AbstractOptimizer, op::HybridProblem)
    This function performs outer approximation algorithm.
"""
function oa_algorithm(model::MOI.AbstractOptimizer, op::OriginalProblem)
    #create linear approximation (milp) model from the original problem
    oad = OAdata()
    init_oa_data(op, oad)
    construct_linear_model(model, op, oad)
    add_oa_cut(model, op, oad)
    #create nlp model from the original problem
    construct_nlp_model(model, op, oad)
    #create reformulated nlp model from the original problem
    construct_ref_nlp_model(model, op, oad)
    reserve_ref_linear_constraints(oad)
    #create linear approximation (milp) model from the reformulated nlp model
    construct_ref_mip_model(model, op, oad)
    #create reformulated feasibility problem
    construct_ref_feasibility_model(oad::OAdata)
    #start of outer approximation algoirthm
    start_time = time()
    oad.oa_started = true
    duration = 10000 #this could be provided as parameter
    while (time() - start_time) < duration
        oad.oa_iter += 1
        #solve reformulated nlp problem
        solve_ref_nlp_model(model, oad);
        oad.nlp_infeasible && solve_ref_feasibility_model(model, oad)
        #update milp master problem
        add_ref_oa_cut(model, oad)
        #solve milp problem
        solve_ref_mip_model(model, oad);
        #check termination conditions
        check_termination_condition(oad)
        if (oad.oa_status == :Optimal) || (oad.oa_status == :Suboptimal) || (oad.oa_status == :FailedOA)
            break
        end
    end
    end_time = time()
    oad.total_time = (end_time-start_time)
    print("\n The outer approximation algorithm has been terminated. \n")
    print("\n The solution status :", oad.oa_status, "\n")
    print("\n The number of oa iterations : ", oad.oa_iter, "\n")
    print("\n The solution time : ", oad.total_time, " secs \n")
    print("\n The objective value :", oad.obj_val, "\n")
    print("\n The objective gap :", oad.obj_gap, "\n")
    print("\n The solution :", oad.incumbent, "\n")
    print("\n The mip solution :", oad.ref_mip_solution, "\n")
end

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
using HerbLanguage: c_const, c_var, c_functor, c_pred, c_type, Variable, Term, List, LPair, Literal, Negation, Constant, Predicate, Structure, Functor, Clause, Conj, get_arity

"""
    Julia -> SWIPL
"""

function to_swipl_ref(c::Constant, previous_variables::Dict{Variable,Cint}, ref::Cint)
    swipl_put_atom_chars(ref, c.name)
end

function to_swipl(c::Constant, previous_variables::Dict{Variable,Cint})
    ref = swipl_new_term_ref()
    to_swipl_ref(c, previous_variables, ref)
    ref
end

function to_swipl_ref(v::Variable, previous_variables::Dict{Variable,Cint}, ref::Cint)
    if haskey(previous_variables, v)
        swipl_put_term(ref, previous_variables[v])
    else
        swipl_put_variable(ref)
        previous_variables[v] = ref
    end
end

function to_swipl(v::Variable, previous_variables::Dict{Variable,Cint})
    if haskey(previous_variables, v)
        previous_variables[v]
    else
        tmpv = swipl_new_term_ref()
        swipl_put_variable(tmpv)
        previous_variables[v] = tmpv
        tmpv
    end
end

function to_swipl_ref(i::Int64, previous_variables::Dict{Variable,Cint}, ref::Cint)
    swipl_put_integer(ref, i)
end

function to_swipl(i::Int64, previous_variables::Dict{Variable,Cint})
    ref = swipl_new_term_ref()
    to_swipl_ref(i, previous_variables, ref)
    ref
end

function to_swipl_ref(i::Float64, previous_variables::Dict{Variable,Cint}, ref::Cint)
    swipl_put_float(ref, i)
end

function to_swipl(i::Float64, previous_variables::Dict{Variable,Cint})
    ref = swipl_new_term_ref()
    to_swipl_ref(i, previous_variables, ref)
    ref
end

function to_swipl(f::Union{Functor,Predicate}, previous_variables::Dict{Variable,Cint})
    func_atm = swipl_new_atom(f.name)
    func = swipl_new_functor(func_atm, f.arity)
    func
end

function to_swipl_ref(s::Union{Structure,Literal}, previous_variables::Dict{Variable, Cint}, ref::Cint)
    if isa(s, Structure)
        func = to_swipl(s.functor, previous_variables)
    else
        func = to_swipl(s.predicate, previous_variables)
    end

    compound_arg = swipl_new_term_refs(get_arity(s))#(s.functor.arity)
    struct_args = s.arguments
    to_swipl_ref(struct_args[begin], previous_variables, compound_arg)
    for ind in 1:(length(struct_args) - 1)
        to_swipl_ref(struct_args[begin+ind], previous_variables, convert(Int32, compound_arg+ind))
    end

    swipl_cons_functor(ref, func, compound_arg)
end

function to_swipl(s::Union{Structure,Literal}, previous_variables::Dict{Variable, Cint})
    structure = swipl_new_term_ref()
    to_swipl_ref(s, previous_variables, structure)
    structure
end

function to_swipl_ref(l::List, previous_variables::Dict{Variable, Cint}, ref::Cint)
    clist_term = swipl_new_term_ref()

    swipl_put_nil(ref)
    for ind in 0:length(l.elements) - 1
        to_swipl_ref(l.elements[end-ind], previous_variables, clist_term)
        swipl_cons_list(ref, clist_term, ref)
    end
end

function to_swipl(l::List, previous_variables::Dict{Variable, Cint})
    list_ref = swipl_new_term_ref()
    to_swipl_ref(l, previous_variables, list_ref)
    list_ref
end

function to_swipl_ref(p::LPair, previous_variables::Dict{Variable,Cint}, ref::Cint)
    head = swipl_new_term_ref()
    tail = swipl_new_term_ref()

    to_swipl_ref(p.head, previous_variables, head)
    to_swipl_ref(p.tail, previous_variables, tail)

    swipl_cons_list(ref, head, tail)
end

function to_swipl(p::LPair, previous_variables::Dict{Variable, Cint})
    ref = swipl_new_term_ref()
    to_swipl_ref(p, previous_variables, ref)
    ref
end

function to_swipl(n::Negation, previous_variables::Dict{Variable, Cint})
    lit_ref = to_swipl(n.literal, previous_variables)

    n_atom = swipl_new_term_ref()
    swipl_put_atom_chars(n_atom, "\\+")

    neg_functor = swipl_new_functor(n_atom, 1)

    neg_ref = swipl_new_term_atom()
    swipl_cons_functor(neg_ref, neg_functor, lit_ref)

    neg_ref
end

function conjoin_swipl(ls::Vector{Cint})
    if length(ls) == 1
        ls[begin]
    else
        c_atom = swipl_new_atom(",")
        conj_functor = swipl_new_functor(c_atom, 2)

        compound_arg = swipl_new_term_refs(2)
        conj = swipl_new_term_ref()
        swipl_put_term(compound_arg, ls[begin])
        swipl_put_term(convert(Int32, compound_arg + 1), conjoin_swipl(ls[begin+1:end]))
        swipl_cons_functor(conj, conj_functor, compound_arg)

        conj
    end
end

function to_swipl(cl::Clause, previous_variables::Dict{Variable,Cint})

    head_lit = to_swipl(cl.head, previous_variables)
    body_lits = [to_swipl(x, previous_variables) for x in cl.body.elements]
    body_terms = conjoin_swipl(body_lits)

    clause_atom = swipl_new_atom(":-")
    clause_functor = swipl_new_functor(clause_atom, 2)

    entire_clause = swipl_new_term_ref()
    compound_args = swipl_new_term_refs(2)
    swipl_put_term(compound_args, head_lit)
    swipl_put_term(convert(Int32, compound_args + 1), body_terms)
    swipl_cons_functor(entire_clause, clause_functor, compound_args)

    entire_clause
end

"""
    SWIPL -> Julia
"""

function swipl_to_const(term::Cint)
    name = swipl_get_atom_chrs(term)
    if islowercase(name[begin])
        c_const(name)
    else
        c_const("\"$(name)\"")
    end
end

function swipl_to_int(term::Cint)
    convert(Int64, swipl_get_integer(term))
end

function swipl_to_float(term::Cint)
    convert(Float64, swipl_get_float(term))
end

function swipl_to_var(term::Cint, swipl_term_to_var::Dict{Cint,Variable})
    if haskey(swipl_term_to_var, term)
        swipl_term_to_var[term]
    else
        existing_names = Set(x.name for x in values(swipl_term_to_var))
        available_names = [x for x in 'A':'Z' if !in(x, existing_names)]
        if length(available_names) == 0
            available_names = ["$x$y" for x in 'A':'Z' for y in 'A':'Z' if !in("$x$y", existing_names)]
        end

        new_name = available_names[begin]
        new_var = c_var(new_name)
        swipl_term_to_var[term] = new_var
        new_var
    end
end

function swipl_to_list(term::Cint)
    elements = Vector{Term}()
    list = swipl_copy_term_ref(term)

    head = swipl_new_term_ref()

    while swipl_get_list(list, head, list)
        push!(from_swipl(head), elements)
    end

    List(elements)
end

function swipl_to_pair(term::Cint)
    head = swipl_new_term_ref()
    tail = swipl_new_term_ref()

    swipl_get_head(term, head)
    swipl_get_tail(term, tail)

    LPair(from_swipl(head), from_swipl(tail))
end

function swipl_to_structure(term::Cint, term_to_var::Dict{Cint,Variable})
    name, arity = swipl_get_name_arity(term)
    functor_name = swipl_atom_chars(name)
    functor = c_functor(functor_name, convert(Int64, arity))

    structure_elements = Vector{Term}()
    for arg_index in 1:arity
        elem = swipl_new_term_ref()
        swipl_get_arg(arg_index, term, elem)
        push!(structure_elements, from_swipl(elem,term_to_var))
    end

    Structure(functor, structure_elements)
end

function from_swipl(term::Cint, term_to_var::Dict{Cint,Variable})
    if swipl_is_atom(term)
        swipl_to_const(term)
    elseif swipl_is_integer(term)
        swipl_to_int(term)
    elseif swipl_is_float(term)
        swipl_to_float(term)
    elseif swipl_is_list(term)
        swipl_to_list(term)
    elseif swipl_is_compound(term)
        swipl_to_structure(term, term_to_var)
    elseif swipl_is_variable(term)
        swipl_to_var(term, term_to_var)
    else
        error("unknown term type $(swipl_term_type(term))")
    end
end

from_swipl(term::Cint) = from_swipl(term, Dict{Cint,Variable}())

"""
    High-level interface
"""

function start(prolog::Swipl)
    swipl_start()
    prolog.initiated = true
end

function stop(prolog::Swipl)
    swipl_stop()
    prolog.initiated = false
end

function cleanup(prolog::Swipl)
     swipl_cleanup()
     prolog.initiated = false
end

is_initialised(prolog::Swipl) = swipl_is_initialised()

function abstract_assert(predicate::String, prolog::Swipl, clause::Union{Clause,Literal})
    vars = Dict{Variable, Cint}()
    swipl_repr = to_swipl(clause, vars)

    asserta_pred = swipl_predicate(predicate, 1)
    query = swipl_open_query(asserta_pred, swipl_repr)
    r = swipl_next_solution(query)
    swipl_close_query(query)

    r
end

asserta(prolog::Swipl, clause::Union{Clause,Literal}) = abstract_assert("asserta", prolog, clause)
assertz(prolog::Swipl, clause::Union{Clause,Literal}) = abstract_assert("assertz", prolog, clause)

function retract(prolog::Swipl, clause::Union{Clause,Literal})
    vars = Dict{Variable, Cint}()
    swipl_repr = to_swipl(clause, vars)

    retract_pred = swipl_predicate("retract", 1)
    query = swipl_opn_query(retract_pred, swipl_repr)
    r = swipl_next_solution(query)
    swipl_close_query(query)

    r
end

function query(prolog::Swipl, literal::Literal; max_solutions::Int=-1)
    query_vars = Dict{Variable,Cint}()
    query_args = swipl_new_term_refs(literal.predicate.arity)

    for ind in 0:literal.predicate.arity-1
        to_swipl_ref(literal.arguments[begin+ind], query_vars, convert(Int32, query_args+ind))
    end

    term_to_var_map = Dict((value, key) for (key, value) in query_vars)

    pred = swipl_predicate(literal.predicate.name, literal.predicate.arity)
    query = swipl_open_query(pred, query_args)

    r = swipl_next_solution(query)

    all_solutions = Vector{Dict{Variable,Union{Term,Int64,Float64}}}()

    while r == 1 & (max_solutions != 0)
        max_solutions -= 1

        current_solution = Dict{Variable,Union{Term,Int64,Float64}}()
        for (var, value) in query_vars
            current_solution[var] = from_swipl(value, term_to_var_map)
        end

        push!(all_solutions, current_solution)

        r = swipl_next_solution(query)
    end

    swipl_close_query(query)
    all_solutions
end

has_solution(prolog::Swipl, literal::Literal) = !isempty(query(prolog, literal; max_solutions=1))


text:
module EmacsREPL

using FileIO
using FileIO: @format_str
using Images: channelview, colorview, RGB
# using Printf

export viewrepl

greet() = print("Hello World!")

# I have to import it before extending it, this is used for debugging
# import EmacsREPL.viewrepl

function viewrepl(img::Array{T,N}) where {T,N}
    if size(img)[3] == 3
        getRGB(X) = colorview(RGB, permutedims(X, (3,1,2)))
        rgbimg = getRGB(img);
        raw_viewrepl(rgbimg)
    elseif size(img)[3] == 1
        raw_viewrepl(img)
    else
        error("Unsupported channel size: ", size(img)[3])
    end
end

# function viewrepl(img)
#     raw_viewrepl(img)
# end

# I need to customize based on the type of img
function raw_viewrepl(img)
    # save("/tmp/a.png")
    path, io = Base.Filesystem.mktemp()
    # FileIO.save(FileIO.Stream(format"PNG", io), img)
    # I'm using .png suffix to hopefully make Emacs' life easier
    # The io variable is ignored, merely used for avoiding race condtion
    FileIO.save(File(format"PNG", path * ".png"), img)
    println("$(path).png")
    println("#<Image: $(path).png>")

    # If saved without .png suffix, and using the path, it is saved as
    # "imagemagick format", and cannot be displayed by 'png, but
    # 'imagemagick will be fast.
    #
    # FileIO.save(File(format"PNG", path), img)

    # If otherwise saved by using io variable, the image is saved in
    # 'png format, and imagemagick will be very slow. So just do NOT
    # use this.
    #
    # save(Stream(format"PNG", io), img)

    # s = @sprintf("#<Image: %s>", path * ".png")

    # The string interpolation just works
    # println(path)
    # println("#<Image: $(path)>")
    # what to return?
    nothing
end

end # module

max_stars_repo_name: lihebi/julia-repl

max_issues_repo_name: lihebi/julia-repl

max_forks_repo_name: lihebi/julia-repl


text:
using Test
using QuantumPropagators
using QuantumControlBase
using LinearAlgebra
using Distributions
using QuantumControlBase.TestUtils

@testset "TLS dissipation" begin

    function ket(i)
        Ψ = zeros(ComplexF64, 2)
        Ψ[i+1] = 1.0
        return Ψ
    end

    bra(i) = adjoint(ket(i))

    ketbra(i, j) = ket(i) * bra(j)

    γ₁ = 0.5
    γ₂ = 0.2

    Âᵧ₁ = √γ₁ * ketbra(0, 1)
    Âᵧ₂ = √(2γ₂) * ketbra(1, 1)

    Ψ₀ = (ket(0) + ket(1)) / √2
    ρ⃗₀ = reshape(Ψ₀ * Ψ₀', :)
    ℒ = Array(liouvillian(nothing, [Âᵧ₁, Âᵧ₂]; convention=:TDSE)[1])

    T = 1.0
    tlist = [0.0, T]

    ρ⃗_out = propagate(ρ⃗₀, (tlist, i; kwargs...) -> ℒ, tlist; method=:expprop)
    ρ_out = reshape(ρ⃗_out, 2, 2)

    #! format: off
    ρ_expected = 0.5 * ComplexF64[
            (2-exp(-γ₁*T))    (exp(-(γ₁/2 + γ₂)*T));
        (exp(-(γ₁/2 + γ₂)*T))       (exp(-γ₁*T))
    ]
    #! format: on

    @test abs(1 - tr(ρ_out)) < 1e-15 # total population
    @test abs(tr(ρ_out^2)) < 1.0  # purity
    @test norm(ρ_out - ρ_expected) < 1e-15

end

@testset "LvN" begin

    N = 100

    function ket(i)
        Ψ = zeros(ComplexF64, N)
        Ψ[i] = 1.0
        return Ψ
    end

    bra(i) = adjoint(ket(i))

    ketbra(i, j) = ket(i) * bra(j)

    Ĥ₀ = random_hermitian_matrix(N, 1)
    Ĥ₁ = random_hermitian_matrix(N, 0.1)

    H = (Ĥ₀, (Ĥ₁, t -> 1.0))
    Ĥ = H[1] + H[2][1] * H[2][2](0)

    Ψ₀ = random_state_vector(N)
    ρ₀ = Ψ₀ * Ψ₀'
    ρ⃗₀ = reshape(ρ₀, :)
    𝕚 = 1im

    ℒ_nodiss = liouvillian(Ĥ; convention=:LvN)[1]
    @test norm(𝕚 * (Ĥ * ρ₀ - ρ₀ * Ĥ) - reshape(ℒ_nodiss * ρ⃗₀, N, N)) < 1e-15

    ℒ_nodiss = liouvillian(Ĥ; convention=:TDSE)[1]
    @test norm((Ĥ * ρ₀ - ρ₀ * Ĥ) - reshape(ℒ_nodiss * ρ⃗₀, N, N)) < 1e-15

    γ₁ = 0.2
    decay_to_ground = [√γ₁ * ketbra(1, i) for i ∈ 2:N]

    γ₂ = 0.1
    dephasing = [√γ₂ * ketbra(i, i) for i ∈ 1:N]

    c_ops = (decay_to_ground..., dephasing...)
    L = liouvillian(H, c_ops; convention=:LvN)
    ℒ = L[1] + L[2][1] * L[2][2](0)

    L = liouvillian(H, c_ops; convention=:LvN)
    ℒ = L[1] + L[2][1] * L[2][2](0)

    ρ̇_LvN = (
        𝕚 * (Ĥ * ρ₀ - ρ₀ * Ĥ) +
        sum([(A * ρ₀ * A' - (A' * A * ρ₀) / 2 - (ρ₀ * A' * A) / 2) for A ∈ c_ops])
    )

    ρ̇ = reshape(ℒ * ρ⃗₀, N, N)

    @test norm(ρ̇ - ρ̇_LvN) < 1e-15

end


text:
#using SurfaceGeometry
using Optim

struct Zinchenko2013{T<:AbstractFloat}
    gamma::AbstractFloat
    C::AbstractFloat
    h2::Vector{T}
    ftol::AbstractFloat
end

function Zinchenko2013(points,faces,n;gamma=0.25,Cp=0.4,ftol=1e-6)

    function VArea(vi)
        area = 0
        for tri in FaceVRing(vi,faces)
            v0,v1,v2 = faces[:,tri]
            y0 = points[:,v0]
            y1 = points[:,v1]
            y2 = points[:,v2]

            Svec = cross(y1-y0,y2-y0)/2
            triarea = norm(Svec)

            area += triarea
        end
        return area
    end

    Lambda = Array{Float64}(undef,size(points,2))
    for i in 1:size(points,2)

        normal = n[:,i]

        p0 = points[:,i]
        vects = Array{Float64}(undef,3,0)
        for vi in VertexVRing(i,faces)
            vects = [vects (points[:,vi]-p0)]
        end

        R = eye(3)
        C,D,E = ZinchenkoDifertential!(vects,normal,R)

        A = [C D/2;D/2 E]
        k1,k2 = eigvals(-A)*2

        lambda = k1^2 + k2^2
        Lambda[i] = lambda
    end

    s = 0
    #gamma = 0.25
    N = size(points,2)
    for i in 1:size(points,2)
        s += Lambda[i]^gamma * VArea(i)/3
    end
    K = 4/sqrt(3)/N * s

    h2 = Array{Float64}(undef,size(points,2))
    for i in 1:size(points,2)
        h2[i] = K*Lambda[i]^(-gamma)
    end

    Zinchenko2013(gamma,Cp,h2,ftol)
end

function F(points,faces,v,zc::Zinchenko2013)

    Cp = zc.C
    h2 = zc.h2

    s = 0
    for ti in 1:size(faces,2)
        v1,v2,v3 = faces[:,ti]

        a = norm(points[:,v2] - points[:,v1])
        b = norm(points[:,v3] - points[:,v2])
        c = norm(points[:,v1] - points[:,v3])

        xava = dot(points[:,v2]-points[:,v1],v[:,v2]-v[:,v1])
        xbvb = dot(points[:,v3]-points[:,v2],v[:,v3]-v[:,v2])
        xcvc = dot(points[:,v1]-points[:,v3],v[:,v1]-v[:,v3])

        s += v2>v1 ? 4*(2/(h2[v2]+h2[v1]) - (h2[v2] + h2[v1])/2/a^4)^2 * xava^2 : 0
        s += v3>v2 ? 4*(2/(h2[v2]+h2[v3]) - (h2[v2] + h2[v3])/2/b^4)^2 * xbvb^2 : 0
        s += v1>v3 ? 4*(2/(h2[v1]+h2[v3]) - (h2[v1] + h2[v3])/2/c^4)^2 * xcvc^2 : 0

        Cdelta = 1/4*sqrt(1 - 2* (a^4 + b^4 + c^4)/(a^2 + b^2 + c^2)^2)
        A = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( a^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        B = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( b^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        C = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( c^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        DCdelta = -A*xava - B*xbvb - C*xcvc

        s += Cp*DCdelta^2/Cdelta^2
    end

    return s
end

function gradF!(points,faces,v,storage,zc::Zinchenko2013)

    Cp = zc.C
    h2 = zc.h2

    storage[:,:] .= 0

    for i in 1:size(faces,2)
        v1,v2,v3 = faces[:,i]
        a = norm(points[:,v2] - points[:,v1])
        b = norm(points[:,v3] - points[:,v2])
        c = norm(points[:,v1] - points[:,v3])

        xava = dot(points[:,v2]-points[:,v1],v[:,v2]-v[:,v1])
        xbvb = dot(points[:,v3]-points[:,v2],v[:,v3]-v[:,v2])
        xcvc = dot(points[:,v1]-points[:,v3],v[:,v1]-v[:,v3])

        Cdelta = 1/4*sqrt(1 - 2 * (a^4 + b^4 + c^4)/(a^2 + b^2 + c^2)^2)
        A = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( a^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        B = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( b^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        C = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( c^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        DCdelta = -A*xava - B*xbvb - C*xcvc

        ha2 = (h2[v2] + h2[v1])/2
        hb2 = (h2[v3] + h2[v2])/2
        hc2 = (h2[v1] + h2[v3])/2

        A_ = 2*Cp*DCdelta*A/Cdelta^2
        B_ = 2*Cp*DCdelta*B/Cdelta^2
        C_ = 2*Cp*DCdelta*C/Cdelta^2

        xa = points[:,v2] - points[:,v1]
        xb = points[:,v3] - points[:,v2]
        xc = points[:,v1] - points[:,v3]

        storage[:,v1] += (A_ - 8*(1/ha2 - ha2/a^4)^2*xava)*xa - C_*xc
        storage[:,v2] += (B_ - 8*(1/hb2 - hb2/b^4)^2*xbvb)*xb - A_*xa
        storage[:,v3] += (C_ - 8*(1/hc2 - hc2/c^4)^2*xcvc)*xc - B_*xb

    end
end

function stabilise!(points,faces,n,v,zc::Zinchenko2013; op=:optim)
    if op==:optim
        stabiliseV2Optim!(points,faces,n,v,zc)
    elseif op===:zinchenko
        stabiliseV2!(points,faces,n,v,zc)
    else
        error("Your requested method is not available")
    end
end

function stabiliseV2Optim!(points,faces,n,v,zc::Zinchenko2013)
    #h2 = hvec(points,faces,n)

    function f(x::Vector)
        vv = reshape(x,3,div(length(x),3))
        return F(points,faces,vv,zc)
    end

    function g!(storage::Vector,x::Vector)

        vv = reshape(x,size(points)...)
        gradf = reshape(storage,size(points)...)

        gradF!(points,faces,vv,gradf,zc)

        for i in 1:size(v,2)
            P = eye(3)-n[:,i]*n[:,i]'
            gradf[:,i] = P*gradf[:,i]
        end
    end

    res = optimize(f,g!,v[:],BFGS())
    #v[:,:] = reshape(res.minimum, size(v)...)[:,:]
    v[:,:] = reshape(Optim.minimizer(res), size(v)...)[:,:]
end

function stabiliseV2!(points,faces,n,v,zc::Zinchenko2013)

    gradF_v = Array{Float64}(undef,size(points)...)
    gradF!(points,faces,v,gradF_v,zc)

    f = Array{Float64}(undef,size(points)...)
    for i in 1:size(points,2)
        f[:,i] = - (eye(3) - n[:,i]*n[:,i]') * gradF_v[:,i]
    end

    Fs = Inf
    eps = 1e-4

    gradF_f = Array{Float64}(undef,size(points)...)
    gradF!(points,faces,f,gradF_f,zc)

    q = 2*dot(gradF_f[:],f[:])
    s = dot(gradF_v[:],f[:]) + dot(gradF_f[:],v[:])
    ksi = - s/q

    gradF_dv = Array{Float64}(undef,size(points)...)
    vp = copy(v)
    v[:,:] = v[:,:] + ksi*f[:,:]

    # step = 1
    for step in 1:150

        gradF!(points,faces,v,gradF_v,zc)
        for i in 1:size(points,2)
            f[:,i] = -(eye(3) - n[:,i]*n[:,i]') * gradF_v[:,i]
        end

        gradF!(points,faces,f,gradF_f,zc)
        gradF!(points,faces,v-vp,gradF_dv,zc)

        S = [dot(gradF_v[:],f[:]) + dot(gradF_f[:],v[:]), dot(gradF_dv[:],v[:]) + dot(gradF_v[:],v[:] - vp[:])]
        Q = [2*dot(gradF_f[:],f[:])   dot(gradF_f[:],v[:]-vp[:])+dot(gradF_dv[:],f[:]); dot(gradF_dv[:],f[:])+dot(gradF_v[:],v[:]-vp[:])   2*dot(gradF_dv[:],v[:]-vp[:])]
        ksi, eta = -Q\S

        v[:,:],vp[:,:] = v[:,:] + ksi*f[:,:]  + eta*(v[:,:] - vp[:,:]), v[:,:]

        Fsnew = F(points,faces,v,zc)
        if (abs(Fs - Fsnew) < eps*Fsnew )
            break
        else
            Fs = Fsnew
        end
    end

    # if step>100
    #     warn("stabilistation used $step steps")
    # end
end


text:
using DashBootstrapComponents

button_group = dbc_buttongroup([
    dbc_button("Left", color = "danger"),
    dbc_button("Middle", color = "warning"),
    dbc_button("Right", color = "success"),
])

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
export sinint, cosint

function sinint{T<:BigFloat}(z::T;n::Integer=2^6)
    Tπ=convert(T,π);dDE=Tπ/2;ga=one(T);b2=one(T);h=log(Tπ*dDE*ga*n/b2)/ga/n
    phi(t::T) = log(exp(Tπ/2*sinh(t))+one(T))
    phip(t::T) = Tπ/2*cosh(t)/(one(T)+exp(-Tπ/2*sinh(t)))
    phiinv(x::T) = asinh(2/Tπ*log(expm1(x)))
    val1,val2=zero(T),zero(T)
    for j=-n:n
        jh=j*h
        phijh = phi(jh)
        temp = exp(-phijh)*phip(jh)/(phijh^2+z^2)*(one(T)-(-1)^(one(T)*j)*exp(-phiinv(abs(z))))
        val1 += temp*z
        val2 += temp*z*phijh
    end
    retval = Tπ/2 - h*(val1*cos(z)+val2*sinc(z/Tπ))
    return z == zero(T) ? zero(T) : z >= zero(T) ? retval : retval-Tπ
end

function Base.sinc{T<:Number}(n::Integer,x::T)
#This program computes the nth sinc differentiation matrices.
    val = zero(T)
    if n ≥ 0
        if x == 0
            val = isodd(n) ? zero(T) : (-1)^(n/2)*convert(T,π)^n/(n+1)
        else
            sp,cp,πk,xk,den = sinpi(x),cospi(x),one(T)/π,one(T),1
            for k=0:n
                m = mod(k,4)
                if m == 0
                    val -= sp*πk*xk/den
                elseif m == 1
                    val += cp*πk*xk/den
                elseif m == 2
                    val += sp*πk*xk/den
                else
                    val -= cp*πk*xk/den
                end
                πk,xk,den = πk*π,xk*x,den*(k+1)
            end
            val *= iseven(n) ? -den/(n+1)/xk : den/(n+1)/xk
        end
    elseif n == -1
        val = 0.5+sinint(π*x)/π
    end
    return val
end
Base.sinc{T<:Number}(n::T,x::T) = sinc(rount(Int,n),x)


text:
export AbstractChain
export ChainConfig
export samples, meta, logq, logw, info, pushsample!
export newchain, step!, sample!
import Base

abstract type AbstractChain{T} <: AbstractVector{T} end

abstract type ChainConfig{T} end

"""
    resizablefields(c::AbstractChain) -> Tuple{Symbol}

Return a tuple of the fieldnames `s` for which `getfield(c, s)` is a vector
which one element for each sample.
"""
function resizablefields(::AbstractChain) end

"""
    function pushsample!(::AbstractChain, sample...) -> AbstractChain

Push a new sample to the chain. This is not `push!` because the sample may be
represented across multiple arguments.
"""
function pushsample!(::AbstractChain, sample...) end

"""
    newchain(config::ChainConfig{T}, args...) -> C

Initialize a new chain of type C.
"""
function newchain(config::ChainConfig, args...) end

"""
    step!(ch::AbstractChain)

Take a single step in the chain, returning the next sample. The return type
depends on AbstractChain subtype. The only change in the chain `ch` is the
iterator state; the samples are unchanged.
"""
function step!(ch::AbstractChain) end

"""
    sample!(ch::AbstractChain, n::Int)

Draw `n` samples from the chain `ch` and append them.
"""
function sample!(ch::AbstractChain, n::Int) end

"""
    getlogdensity(ch::AbstractChain{T})

Get the log-density function the samples are drawn from. Returns a `T -> Real`
function.
"""
function getlogdensity(chain::AbstractChain) end

using MappedArrays

meta(ch::AbstractChain) = getfield(ch, :meta)
globals(ch::AbstractChain) = getfield(ch, :globals)
logw(chain::AbstractChain) = mappedarray(_ -> 0.0, logq(chain))
info(chain::AbstractChain) = getfield(chain, :info)

"""
    samples(chain::AbstractChain{T}) -> AbstractVector{T}

Return a vector of the samples in the chain.
"""
samples(chain::AbstractChain) = getfield(chain, :samples)

"""
    logq(chain::AbstractChain) -> AbstractVector{AbstractFloat}

Return a vector of the log-densities of the sampled chain. Should satisfy

    logq(chain) == getlogdensity(chain).(samples(chain))
"""
logq(chain::AbstractChain) = getfield(chain, :logq)

Base.size(ch::AbstractChain) = size(samples(ch))

Base.length(ch::AbstractChain) = length(samples(ch))

Base.getindex(ch::AbstractChain, n) = getindex(samples(ch), n)

Base.propertynames(ch::AbstractChain) = propertynames(first(samples(ch)))

# TODO: Make this pass @code_warntype
Base.getproperty(ch::AbstractChain, k::Symbol) = getproperty(samples(ch), k)

function Base.showarg(io::IO, chain::AbstractChain{T}, toplevel) where T
    io = IOContext(io, :compact => true)
    print(io, "Chain")
    toplevel && println(io, " with schema ", schema(T))
end

function Base.show(io::IO, ::MIME"text/plain", chain::AbstractChain)
    summary(io, chain)
    print(io, summarize(chain))
end

function Base.show(io::IO, ::MIME"text/html", chain::AbstractChain)
    summary(io, chain)
    print(io, summarize(chain))
end

function Base.resize!(ch::AbstractChain, n::Int)
    for v in resizablefields(ch)
        resize!(v, n)
    end
end

NestedTuples.rmap(f, ch::AbstractChain) = rmap(f, samples(ch))


text:
const MONITORING = Ref{Bool}(false)
const MONITORING_TASK = Ref{Union{Nothing,Task}}(nothing)

"Checks if we are currently monitoring."
ismonitoring() = MONITORING[]
_monitoring!(val::Bool) = MONITORING[] = val
_set_monitoring_task(t::Task) = MONITORING_TASK[] = t
_get_monitoring_task() = MONITORING_TASK[]

struct MonitoringResults
    times::Vector{Float64}
    devices::Vector{Tuple{String,UUID}}
    results::Dict{Symbol,Vector{Vector{Float64}}}
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, r::MonitoringResults)
    println(io, "MonitoringResults")
    return print(io, "└ Quantities: ", collect(keys(r.results)))
end
function Base.show(io::IO, r::MonitoringResults)
    return print(io, "MonitoringResults(", collect(keys(r.results)), ")")
end

function Base.getproperty(r::MonitoringResults, s::Symbol)
    if hasfield(typeof(r), s)
        return getfield(r, s)
    else
        return r.results[s]
    end
end

function Base.:(==)(r1::MonitoringResults, r2::MonitoringResults)
    equal = true
    for field in fieldnames(MonitoringResults)
        if getfield(r1, field) != getfield(r2, field)
            equal = false
        end
    end
    return equal
end

"""
    monitoring_start(; devices=CUDA.devices(), verbose=true)
Start monitoring of GPU temperature, utilization, power usage, etc.

**Keyword arguments:**
* `freq` (default: `1`): polling rate in Hz.
* `devices` (default: `CUDA.devices()`): `CuDevice`s or `NVML.Device`s to monitor.
* `thread` (default: `Threads.nthreads()`): id of the Julia thread that should run the monitoring.
* `verbose` (default: `true`): toggle verbose output.

See also [`monitoring_stop`](@ref).
"""
function monitoring_start(;
    freq=1, devices=CUDA.devices(), thread=Threads.nthreads(), verbose=true
)
    if ismonitoring()
        error("We are already monitoring.")
    end
    if thread <= 0 || thread > Threads.nthreads()
        throw(
            ArgumentError(
                "Invalid value for `thread` (must be `1 ≤ thread ≤ Threads.nthreads()`)."
            ),
        )
    end
    if Threads.nthreads() < length(devices) + 1
        @warn(
            "Notice that there aren't enough Julia threads to test all available GPUs and monitoring them at the same time."
        )
    end
    if eltype(devices) == CuDevice
        devs = _cudevice2nvmldevice.(devices)
    else
        devs = devices
    end
    verbose && @info("Spawning monitoring on Julia thread $thread.")
    # check which quantities we can monitor
    func_symbol_vec = Tuple{Function,Symbol}[]
    if all(supports_get_temperature(dev) for dev in devs)
        push!(func_symbol_vec, (get_temperatures, :temperature))
    end
    if all(supports_get_power_usage(dev) for dev in devs)
        push!(func_symbol_vec, (get_power_usages, :power))
    end
    if all(supports_get_gpu_utilization(dev) for dev in devs)
        push!(
            func_symbol_vec,
            (devs -> getproperty.(get_gpu_utilizations(devs), :compute), :compute),
        )
        push!(
            func_symbol_vec, (devs -> getproperty.(get_gpu_utilizations(devs), :mem), :mem)
        )
    end

    # spawn monitoring
    _monitoring!(true)
    t = @tspawnat thread _monitor(func_symbol_vec; freq, devices=devs)
    _set_monitoring_task(t)
    return nothing
end

function _monitor(func_symbol_vec; freq, devices)
    nfuncs = length(func_symbol_vec)
    t_freq = 1.0 / freq
    # init result dict
    d = Dict{Symbol,Vector{Vector{Float64}}}()
    for i in 1:nfuncs
        f, s = func_symbol_vec[i]
        d[s] = Vector{Float64}[f(devices)]
    end
    times = [0.0]
    # monitor
    t = time()
    while ismonitoring()
        sleep(t_freq)
        push!(times, time() - t)
        for i in 1:nfuncs
            f, s = func_symbol_vec[i]
            push!(d[s], f(devices))
        end
    end
    device_info = Vector{Tuple{String,UUID}}(undef, length(devices))
    for (i, nvmldev) in enumerate(devices)
        dev = _nvmldevice2cudevice(nvmldev)
        device_info[i] = (_device2string(dev), uuid(dev))
    end
    return MonitoringResults(times, device_info, d)
end

"""
    monitoring_stop(; verbose=true) -> results
Stops the GPU monitoring and returns the measured values.
Specifically, `results` is a named tuple with the following keys:
* `time`: the (relative) times at which we measured
* `temperature`, `power`, `compute`, `mem`

See also [`monitoring_start`](@ref) and [`plot_monitoring_results`](@ref).
"""
function monitoring_stop(; verbose=true)::MonitoringResults
    if ismonitoring()
        verbose && @info("Stopping monitoring and fetching results...")
        _monitoring!(false)
        results = fetch(_get_monitoring_task())
        return results
    else
        error("We aren't monitoring, so can't stop it...")
    end
end

"""
    plot_monitoring_results(r::MonitoringResults, symbols=keys(r.results))
Plot the quantities specified through `symbols` of a `MonitoringResults` object.
Will generate a textual in-terminal / in-logfile plot using UnicodePlots.jl.
"""
function plot_monitoring_results(r::MonitoringResults, symbols=keys(r.results))
    for s in symbols
        display(plot_monitoring_results(r, s))
    end
    return nothing
end

function plot_monitoring_results(r::MonitoringResults, s::Symbol)
    println() # top margin
    times = r.times
    values = r.results[s]
    title, ylabel = _symbol2title_and_label(s)
    ylims = _defaultylims(values)
    device_labels = [str for (str, uuid) in r.devices]

    p = UnicodePlots.lineplot(
        times,
        getindex.(values, 1);
        title=title,
        xlabel="Time [s]",
        ylabel=ylabel,
        name=device_labels[1],
        ylim=ylims,
    )
    for i in 2:length(first(values))
        UnicodePlots.lineplot!(p, times, getindex.(values, i); name=device_labels[i])
    end
    return p
end

"""
    livemonitor_temperature(duration) -> times, temperatures

Monitor the temperature of GPU(s) over a given time period, as specified by `duration` (in seconds).
Returns the (relative) times as a `Vector{Float64}` and the temperatures as a `Vector{Vector{Float64}}`.

For general keyword arguments, see [`livemonitor_something`](@ref).
"""
function livemonitor_temperature(duration; kwargs...)
    return livemonitor_something(
        get_temperatures, duration; title="GPU Temperatures", ylabel="T [C]", kwargs...
    )
end

"""
    livemonitor_powerusage(duration) -> times, powerusage

Monitor the power usage of GPU(s) (in Watts) over a given time period, as specified by `duration` (in seconds).
Returns the (relative) times as a `Vector{Float64}` and the power usage as a `Vector{Vector{Float64}}`.

For general keyword arguments, see [`livemonitor_something`](@ref).
"""
function livemonitor_powerusage(duration; kwargs...)
    return livemonitor_something(
        get_power_usages, duration; title="GPU Power Usage", ylabel="P [W]", kwargs...
    )
end

"""
    livemonitor_something(f, duration) -> times, values

Monitor some property of GPU(s), as specified through the function `f`, over a given time period, as specified by `duration` (in seconds).
Returns the (relative) times as a `Vector{Float64}` and the temperatures as a `Vector{Vector{Float64}}`.

The function `f` will be called on a vector of devices (`CuDevice`s or `NVML.Device`s) and should return a vector of `Float64` values.

**Keyword arguments:**
* `freq` (default: `1`): polling rate in Hz.
* `devices` (default: `NVML.devices()`): `CuDevice`s or `NVML.Device`s to consider.
* `plot` (default: `false`): Create a unicode plot after the monitoring.
* `liveplot` (default: `false`): Create and update a unicode plot during the monitoring. Use optional `ylims` to specify fixed y limits.
* `title` (default: `""`): Title used in unicode plots.
* `ylabel` (default: `"Values"`): y label used in unicode plots.

See: [`livemonitor_temperature`](@ref), [`livemonitor_powerusage`](@ref)
"""
function livemonitor_something(
    f::F,
    duration;
    freq=1,
    devices=NVML.devices(),
    plot=false,
    liveplot=true,
    ylims=nothing,
    title="",
    ylabel="Values",
) where {F}
    t_freq = 1.0 / freq
    values = Vector{Float64}[f(devices)]
    times = Float64[0.0]

    t = time()
    while (time() - t) < duration
        sleep(t_freq)
        push!(times, time() - t)
        push!(values, f(devices))
        if liveplot
            print("\033c")
            display(
                _plot_something(
                    values,
                    times;
                    devices,
                    ylims=isnothing(ylims) ? _defaultylims(values) : ylims,
                    title,
                    ylabel,
                ),
            )
        end
    end

    plot && display(
        _plot_something(
            values,
            times;
            devices,
            ylims=isnothing(ylims) ? _defaultylims(values) : ylims,
            title,
            ylabel,
        ),
    )
    return times, values
end

function _plot_something(
    values,
    times=1:length(values);
    devices=nothing,
    ylims=_defaultylims(values),
    title="",
    ylabel="Values",
)
    if isnothing(devices)
        device_labels = ["" for _ in 0:(length(first(values)) - 1)]
    else
        device_labels = [_device2string(dev) for dev in devices]
    end

    p = lineplot(
        times,
        getindex.(values, 1);
        title=title,
        xlabel="Time [s]",
        ylabel=ylabel,
        name=device_labels[1],
        ylim=ylims,
    )
    for i in 2:length(first(values))
        lineplot!(p, times, getindex.(values, i); name=device_labels[i])
    end
    return p
end

function _defaultylims(values)
    total_max = maximum(maximum(vs) for vs in values)
    total_min = minimum(minimum(vs) for vs in values)
    return (floor(Int, total_min * 0.95), ceil(Int, total_max * 1.05))
end

function _symbol2title_and_label(s::Symbol)
    if s == :temperature
        return "GPU Temperature", "T [C]"
    elseif s == :power
        return "GPU Power Usage", "P [W]"
    elseif s == :compute
        return "GPU Utilization (Compute)", "U [%]"
    elseif s == :mem
        return "GPU Utilization (Memory)", "U [%]"
    else
        return "", "Values"
    end
end

"""
    savefig_monitoring_results(r::MonitoringResults, symbols=keys(r.results))
Save plots of the quantities specified through `symbols` of a `MonitoringResults` object to disk.
**Note:** Only available if CairoMakie.jl is loaded next to GPUInspector.jl.
"""
function savefig_monitoring_results(r::MonitoringResults, symbols=keys(r.results); ext=:pdf)
    return error("You need to load CairoMakie.jl first.")
end