text:
@testset "MSAEditing.jl" begin

    msa_types = (
        Matrix{Residue},
        NamedResidueMatrix{Array{Residue,2}},
        MultipleSequenceAlignment,
        AnnotatedMultipleSequenceAlignment
        )

    pf09645_sto = joinpath(pwd(), "data", "PF09645_full.stockholm")
    gaoetal2011 = joinpath(pwd(), "data", "Gaoetal2011.fasta")

    gaoetal_msas = [ read(gaoetal2011, FASTA, T) for T in msa_types ]
    pfam_msas    = [ read(pf09645_sto, Stockholm, T) for T in msa_types ]

    pfam    = pfam_msas[end]
    pfam_na = pfam_msas[end-1] # na: not annotated

    @testset "Boolean mask array" begin

        matrix_mask = pfam[1:1,:] .== GAP
        @test size(matrix_mask) == (1,110)

        @test size(MSA._column_mask(matrix_mask, pfam)) == (110,)
        @test MSA._column_mask(vec(matrix_mask), pfam) == vec(matrix_mask)
        @test MSA._column_mask(col -> col[1] == GAP, pfam) ==
              MSA._column_mask(matrix_mask, pfam)

        matrix_mask = pfam[:,1:1] .== Residue('-')
        @test size(matrix_mask) == (4,1)

        @test size(MSA._sequence_mask(matrix_mask, pfam)) == (4,)
        @test MSA._sequence_mask(vec(matrix_mask), pfam) == vec(matrix_mask)
        @test MSA._sequence_mask(seq -> seq[1] == GAP, pfam) ==
              MSA._sequence_mask(matrix_mask, pfam)
    end

    @testset "filtersequences!" begin

        for msa in pfam_msas[3:4]
            @test filtersequences!(copy(msa), 1:4 .< 3) == filtersequences(msa, 1:4 .< 3)
            @test nsequences(msa) == 4
        end

        for msa in pfam_msas
            @test getsequence(filtersequences(
                msa, [1,2,3,4] .== 3), 1) == getsequence(msa, 3)
            @test getsequence(filtersequences(
                msa, [1,2,3,4] .> 2), 2) == getsequence(msa, 4)

            @test getsequence(filtersequences(
                msa, Bool[false,false,true,true] ), 2) == getsequence(msa, 4)

            @test_throws AssertionError filtersequences(msa, [1,2,3,4,5,6,7,8,9,10] .> 2)
            @test_throws AssertionError filtersequences(msa, [1,2,3] .> 2)
        end
    end

    @testset "filtercolumns!" begin

        for msa in pfam_msas[3:4]
            @test filtercolumns!(copy(msa), 1:110 .< 11) == filtercolumns(msa, 1:110 .< 11)
            @test ncolumns(msa) == 110
        end

        for msa in pfam_msas
            @test vec(getresidues(getsequence(filtercolumns(
                msa, collect(1:110) .<= 10), 4))) == res"QTLNSYKMAE"

            @test_throws AssertionError filtercolumns(msa, [1,2,3] .> 2)
            @test_throws AssertionError filtercolumns(msa, collect(1:200) .<= 10)
        end

        @testset "filtercolumns! for sequences" begin

            annseq = getsequence(pfam, 4)
            seq = getsequence(pfam_na, 4) # na: not annotated

            # Sequences are matrices
            @test annseq[1, vec(annseq .== Residue('Q'))] == res"QQQQQQQQQQQQQQ"
            @test seq[1, vec(seq .== Residue('Q'))] == res"QQQQQQQQQQQQQQ"

            filtered_annseq = filtercolumns!(copy(annseq), annseq .== Residue('Q'))
            filtered_seq = filtercolumns!(copy(seq), seq .== Residue('Q'))

            @test filtercolumns(seq, seq .== Residue('Q')) ==
                filtercolumns(seq, seq .== Residue('Q'))

            @test_throws AssertionError filtercolumns(annseq, 1:(length(seq)-10) .> 2)
            @test_throws AssertionError filtercolumns(seq, 1:(length(seq)+10) .<= 10)

            # Sequences are matrices
            @test vec(getresidues(filtered_annseq)) == res"QQQQQQQQQQQQQQ"
            @test vec(getresidues(filtered_seq)) == res"QQQQQQQQQQQQQQ"

            @test getannotcolumn(filtered_annseq, "SS_cons") == "XHHEXXXXXXXXXX"
            @test getannotresidue(filtered_annseq,
                "F112_SSV1/3-112", "SS") == "XHHEXXXXXXXXXX"
        end
    end

    @testset "Reference" begin

        @testset "setreference" begin

            for msa in pfam_msas[2:4]
                copy_msa = copy(msa)
                setreference!(copy_msa, 4)
                # Sequences are matrices
                @test vec(getresidues(getsequence(copy_msa,1))) ==
                    res"QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ"
                setreference!(copy_msa, "C3N734_SULIY/1-95")
                @test_throws ErrorException setreference!(copy_msa, "FALSE_SEQID/1-95")
                @test vec(getresidues(getsequence(copy_msa,4))) ==
                    res"QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ"
                @test copy_msa == msa
            end
        end

        @testset "gapstrip!" begin

            for msa in pfam_msas[3:4]
                copy_msa = copy(msa)

                setreference!(copy_msa, 4)
                gapstrip!(copy_msa, gaplimit=1.0, coveragelimit=0.0)
                @test size(copy_msa) == (4, 110)

                setreference!(copy_msa, 1)
                gapstrip!(copy_msa)
                residuefraction(copy_msa[1,:]) == 1.0
            end

            for msa in pfam_msas[1:2]
                @test size(gapstrip(msa, gaplimit=1.0, coveragelimit=0.0)) == (4, 92)
                @test residuefraction(gapstrip(msa)[1,:]) == 1.0
            end
        end

        @testset "adjustreference!" begin

            for msa in pfam_msas[3:4]
                copy_msa = copy(msa)
                adjustreference!(copy_msa)
                residuefraction(copy_msa[1,:]) == 1.0
            end

            for msa in pfam_msas[1:2]
                @test residuefraction(adjustreference(msa)[1,:]) == 1.0
            end
        end
    end
end


text:
module ExampleFunctions

using ..DZOptimization: unsafe_sqrt

function rosenbrock_objective(v::AbstractVector{T}) where {T}
    @inbounds x, y = v[1], v[2]
    t1 = 1 - x
    t2 = y - x * x
    return t1 * t1 + 100 * (t2 * t2)
end

function rosenbrock_gradient!(g::AbstractVector{T},
                              v::AbstractVector{T}) where {T}
    @inbounds x, y = v[1], v[2]
    t1 = 1 - x
    t2 = y - x * x
    @inbounds g[1] = -2 * t1 - 400 * x * t2
    @inbounds g[2] = 200 * t2
    return g
end

function riesz_energy(points::AbstractMatrix{T}) where {T}
    dim, num_points = size(points)
    energy = zero(T)
    @inbounds for i = 2 : num_points
        for j = 1 : i-1
            dist_sq = zero(T)
            @simd ivdep for k = 1 : dim
                dist = points[k,i] - points[k,j]
                dist_sq += dist * dist
            end
            energy += inv(unsafe_sqrt(dist_sq))
        end
    end
    return energy
end

function riesz_gradient!(
        grad::AbstractMatrix{T}, points::AbstractMatrix{T}) where {T}
    dim, num_points = size(points)
    @inbounds for j = 1 : num_points
        @simd ivdep for k = 1 : dim
            grad[k,j] = zero(T)
        end
        for i = 1 : num_points
            if i != j
                dist_sq = zero(T)
                @simd ivdep for k = 1 : dim
                    dist = points[k,i] - points[k,j]
                    dist_sq += dist * dist
                end
                inv_dist_cubed = unsafe_sqrt(dist_sq) / (dist_sq * dist_sq)
                @simd ivdep for k = 1 : dim
                    grad[k,j] += (points[k,i] - points[k,j]) * inv_dist_cubed
                end
            end
        end
    end
    return grad
end

function riesz_hessian!(
        hess::AbstractArray{T,4}, points::AbstractMatrix{T}) where {T}
    dim, num_points = size(points)
    @inbounds for k = 1 : num_points
        for s = 1 : num_points
            if s == k
                for l = 1 : dim
                    @simd ivdep for t = 1 : dim
                        hess[t,s,l,k] = zero(T)
                    end
                end
                for j = 1 : num_points
                    if j != k
                        dist_sq = zero(T)
                        @simd ivdep for d = 1 : dim
                            temp = points[d,k] - points[d,j]
                            dist_sq += temp * temp
                        end
                        dist = unsafe_sqrt(dist_sq)
                        dist_cb = dist * dist_sq
                        for l = 1 : dim
                            @simd ivdep for t = 1 : dim
                                plkj = points[l,k] - points[l,j]
                                ptkj = points[t,k] - points[t,j]
                                temp = (plkj * ptkj) / (dist_sq * dist_cb)
                                hess[t,s,l,k] += (temp + temp + temp)
                            end
                            hess[l,s,l,k] -= inv(dist_cb)
                        end
                    end
                end
            else
                dist_sq = zero(T)
                @simd ivdep for d = 1 : dim
                    temp = points[d,k] - points[d,s]
                    dist_sq += temp * temp
                end
                dist = unsafe_sqrt(dist_sq)
                dist_cb = dist * dist_sq
                for l = 1 : dim
                    @simd ivdep for t = 1 : dim
                        plks = points[l,k] - points[l,s]
                        ptsk = points[t,s] - points[t,k]
                        temp = (plks * ptsk) / (dist_sq * dist_cb)
                        hess[t,s,l,k] = (temp + temp + temp)
                    end
                    hess[l,s,l,k] += inv(dist_cb)
                end
            end
        end
    end
    return hess
end

function constrain_riesz_gradient_sphere!(
        grad::AbstractMatrix{T}, points::AbstractMatrix{T}) where {T}
    dim, num_points = size(points)
    @inbounds for i = 1 : num_points
        overlap = zero(T)
        @simd ivdep for k = 1 : dim
            overlap += points[k,i] * grad[k,i]
        end
        @simd ivdep for k = 1 : dim
            grad[k,i] -= overlap * points[k,i]
        end
    end
    return grad
end

function constrain_riesz_hessian_sphere!(
        hess::AbstractArray{T,4}, points::AbstractMatrix{T},
        unconstrained_grad::AbstractMatrix{T}) where {T}
    dim, num_points = size(points)
    @inbounds for s = 1 : num_points
        for t = 1 : num_points
            for i = 1 : dim
                temp = zero(T)
                @simd ivdep for j = 1 : dim
                    temp += hess[j,s,i,t] * points[j,s]
                end
                @simd ivdep for j = 1 : dim
                    hess[j,s,i,t] -= temp * points[j,s]
                end
            end
        end
        for t = 1 : num_points
            for i = 1 : dim
                temp = zero(T)
                @simd ivdep for j = 1 : dim
                    temp += hess[i,t,j,s] * points[j,s]
                end
                @simd ivdep for j = 1 : dim
                    hess[i,t,j,s] -= temp * points[j,s]
                end
            end
        end
    end
    @inbounds for s = 1 : num_points
        for i = 1 : dim
            gis = unconstrained_grad[i,s]
            pis = points[i,s]
            for j = 1 : dim
                pjs = points[j,s]
                for k = 1 : dim
                    pks = points[k,s]
                    temp = pis * pjs * pks
                    temp = temp + temp + temp
                    if j == k; temp -= pis; end
                    if i == k; temp -= pjs; end
                    if i == j; temp -= pks; end
                    hess[k,s,j,s] += gis * temp
                end
            end
        end
    end
    return hess
end

end # module ExampleFunctions

ext: jl


text:
struct Screen
  window
  renderer
  height::Int
  width::Int
  background::ARGB

  function Screen(name, w, h, color)
    win, renderer = make_win_renderer(name, w, h)
    if !(color isa ARGB)
      color = ARGB(color)
    end
    new(win, renderer, h, w, color)
  end
end

function clear(s::Screen)
  fill(s, s.background)
end

clear() = clear(game[].screen)

Base.fill(c::Colorant) = Base.fill(game[].screen, c)

function Base.fill(s::Screen, c::Colorant)
  SDL2.SetRenderDrawColor(s.renderer, sdl_colors(c)...)
  SDL2.RenderClear(s.renderer)
end


text:
function calc_CRB(s, w)
    s = reshape(s, size(s,1)*size(s,2), size(s,3))
    real.(w * diag(inv(s' * s)))
end

function grad_ω1_fd(w, ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list)
    Δω1 = 1e-2
    _grad_ω1 = similar(ω1)

    α = ω1 .* TRF
    s0 = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list=grad_list)
    CRB0 = calc_CRB(s0, w)

    for t in eachindex(ω1)
        α[t] = (ω1[t] + Δω1) * TRF[t]

        ds = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list=grad_list)
        CRB1 = calc_CRB(ds, w)

        α[t] = ω1[t] .* TRF[t]

        _grad_ω1[t] = (CRB1 - CRB0) / Δω1
    end
    _grad_ω1[1] = 0
    return _grad_ω1
end

function grad_TRF_fd(w, ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list)
    ΔTRF = 1e-9
    _grad_TRF = similar(ω1)

    α = ω1 .* TRF
    s0 = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list=grad_list)
    CRB0 = calc_CRB(s0, w)

    for t in eachindex(ω1)
        TRF[t] += ΔTRF
        α[t] = ω1[t] * TRF[t]

        ds = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list=grad_list)
        CRB1 = calc_CRB(ds, w)
        _grad_TRF[t] = (CRB1 - CRB0) / ΔTRF

        TRF[t] -= ΔTRF
        α[t] = ω1[t] * TRF[t]
    end
    _grad_TRF[1] = 0
    return _grad_TRF
end

function dCRBdm_fd(m, w)
    _dCRBdx = similar(m, Float64)
    _dCRBdy = similar(m, Float64)
    Δm = 1e-9

    CRB0 = calc_CRB(m, w)
    for i in eachindex(m)
        dm = copy(m)
        dm[i] = m[i] .+ Δm
        _dCRBdx[i] = real(CRB0 - calc_CRB(dm, w)) / Δm

        dm = copy(m)
        dm[i] = m[i] + 1im * Δm
        _dCRBdy[i] += real(CRB0 - calc_CRB(dm, w)) / Δm
    end

    d(t, r, g) = @SVector [_dCRBdx[t,r,1], _dCRBdy[t,r,1],0,0,0,_dCRBdx[t,r,g + 1], _dCRBdy[t,r,g + 1],0,0,0,0]

    return (CRB0, d)
end

###########################
function grad_TV_ω1_fd(weights, ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, λ_TV)
    Δω1 = 1e-2
    _grad_ω1_fd = similar(ω1)

    (F0, _grad_TV_ω1, _) = MRIgeneralizedBloch.OCT_TV_gradient(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, weights, λ_TV)

    for t in eachindex(ω1)
        ω1[t] += Δω1
        (F1, _, _) = MRIgeneralizedBloch.OCT_TV_gradient(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, weights, λ_TV)
        ω1[t] -= Δω1

        _grad_ω1_fd[t] = (F1 - F0) / Δω1
    end
    _grad_ω1_fd[1] = 0
    return (_grad_TV_ω1, _grad_ω1_fd)
end


text:
# Continuous linear problems

# Basic solver, query, resolve
function linear1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # simple 2 variable, 1 constraint problem
    # min -x
    # st   x + y <= 1   (x + y - 1 ∈ Nonpositives)
    #       x, y >= 0   (x, y ∈ Nonnegatives)

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}),(MOI.SingleVariable,MOI.GreaterThan{Float64})])

    #@test MOI.get(model, MOI.SupportsAddConstraintAfterSolve())
    #@test MOI.get(model, MOI.SupportsAddVariableAfterSolve())
    #@test MOI.get(model, MOI.SupportsDeleteConstraint())

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    v = MOI.addvariables!(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    cf = MOI.ScalarAffineFunction(v, [1.0,1.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf), MOI.LessThan{Float64})
    c = MOI.addconstraint!(model, cf, MOI.LessThan(1.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    vc1 = MOI.addconstraint!(model, MOI.SingleVariable(v[1]), MOI.GreaterThan(0.0))
    # test fallback
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    vc2 = MOI.addconstraint!(model, v[2], MOI.GreaterThan(0.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 2

    objf = MOI.ScalarAffineFunction(v, [-1.0,0.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MinSense

    if config.query
        @test MOI.canget(model, MOI.ListOfVariableIndices())
        vrs = MOI.get(model, MOI.ListOfVariableIndices())
        @test vrs == v || vrs == reverse(v)

        @test !MOI.canget(model, MOI.ObjectiveFunction{MOI.SingleVariable}())
        @test MOI.canget(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test objf ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())

        @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c))
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)

        @test MOI.canget(model, MOI.ConstraintSet(), typeof(c))
        s = MOI.get(model, MOI.ConstraintSet(), c)
        @test s == MOI.LessThan(1.0)

        @test MOI.canget(model, MOI.ConstraintSet(), typeof(vc1))
        s = MOI.get(model, MOI.ConstraintSet(), vc1)
        @test s == MOI.GreaterThan(0.0)

        @test MOI.canget(model, MOI.ConstraintSet(), typeof(vc2))
        s = MOI.get(model, MOI.ConstraintSet(), vc2)
        @test s == MOI.GreaterThan(0.0)
    end

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1, 0] atol=atol rtol=rtol

        @test MOI.canget(model, MOI.ConstraintPrimal(), typeof(c))
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1 atol=atol rtol=rtol

        if config.duals
            @test MOI.canget(model, MOI.DualStatus())
            @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -1 atol=atol rtol=rtol

            # reduced costs
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc1))
            @test MOI.get(model, MOI.ConstraintDual(), vc1) ≈ 0 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc2))
            @test MOI.get(model, MOI.ConstraintDual(), vc2) ≈ 1 atol=atol rtol=rtol
        end
    end

    # change objective to Max +x

    objf = MOI.ScalarAffineFunction(v, [1.0,0.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.query
        @test MOI.canget(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test objf ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    end

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1, 0] atol=atol rtol=rtol

        if config.duals
            @test MOI.canget(model, MOI.DualStatus())
            @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -1 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc1))
            @test MOI.get(model, MOI.ConstraintDual(), vc1) ≈ 0 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc2))
            @test MOI.get(model, MOI.ConstraintDual(), vc2) ≈ 1 atol=atol rtol=rtol
        end
    end

    # add new variable to get :
    # max x + 2z
    # s.t. x + y + z <= 1
    # x,y,z >= 0

    MOI.canaddvariable(model)
    z = MOI.addvariable!(model)
    push!(v, z)
    @test v[3] == z

    if config.query
        # Test that the modifcation of v has not affected the model
        @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c))
        vars = MOI.get(model, MOI.ConstraintFunction(), c).variables
        @test vars == [v[1], v[2]] || vars == [v[2], v[1]]
        @test MOI.canget(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test MOI.ScalarAffineFunction([v[1]], [1.0], 0.0) ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    end

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    vc3 = MOI.addconstraint!(model, MOI.SingleVariable(v[3]), MOI.GreaterThan(0.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 3

    @test MOI.canmodifyconstraint(model, c, MOI.ScalarCoefficientChange{Float64})
    MOI.modifyconstraint!(model, c, MOI.ScalarCoefficientChange{Float64}(z, 1.0))

    @test MOI.canmodifyobjective(model, MOI.ScalarCoefficientChange{Float64})
    MOI.modifyobjective!(model, MOI.ScalarCoefficientChange{Float64}(z, 2.0))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 3

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus(1))
        @test MOI.get(model, MOI.PrimalStatus(1)) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0, 0, 1] atol=atol rtol=rtol

        @test MOI.canget(model, MOI.ConstraintPrimal(), typeof(c))
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1 atol=atol rtol=rtol

        if config.duals
            @test MOI.canget(model, MOI.DualStatus())
            @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -2 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc1))
            @test MOI.get(model, MOI.ConstraintDual(), vc1) ≈ 1 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc2))
            @test MOI.get(model, MOI.ConstraintDual(), vc2) ≈ 2 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc3))
            @test MOI.get(model, MOI.ConstraintDual(), vc3) ≈ 0 atol=atol rtol=rtol
        end
    end

    # setting lb of x to -1 to get :
    # max x + 2z
    # s.t. x + y + z <= 1
    # x >= -1
    # y,z >= 0
    @test MOI.canmodifyconstraint(model, vc1, MOI.GreaterThan{Float64})
    MOI.modifyconstraint!(model, vc1, MOI.GreaterThan(-1.0))

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3 atol=atol rtol=rtol
    end

    # put lb of x back to 0 and fix z to zero to get :
    # max x + 2z
    # s.t. x + y + z <= 1
    # x, y >= 0, z = 0
    @test MOI.canmodifyconstraint(model, vc1, MOI.GreaterThan{Float64})
    MOI.modifyconstraint!(model, vc1, MOI.GreaterThan(0.0))

    @test MOI.candelete(model, vc3)
    MOI.delete!(model, vc3)

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.EqualTo{Float64})
    vc3 = MOI.addconstraint!(model, MOI.SingleVariable(v[3]), MOI.EqualTo(0.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 2

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1 atol=atol rtol=rtol
    end

    # modify affine linear constraint set to be == 2 to get :
    # max x + 2z
    # s.t. x + y + z == 2
    # x,y >= 0, z = 0
    @test MOI.candelete(model, c)
    MOI.delete!(model, c)
    cf = MOI.ScalarAffineFunction(v, [1.0,1.0,1.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf), MOI.EqualTo{Float64})
    c = MOI.addconstraint!(model, cf, MOI.EqualTo(2.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 0
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}}()) == 1

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2 atol=atol rtol=rtol
    end

    # modify objective function to x + 2y to get :
    # max x + 2y
    # s.t. x + y + z == 2
    # x,y >= 0, z = 0

    objf = MOI.ScalarAffineFunction(v, [1.0,2.0,0.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 4 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0, 2, 0] atol=atol rtol=rtol
    end

    # add constraint x - y >= 0 to get :
    # max x+2y
    # s.t. x + y + z == 2
    # x - y >= 0
    # x,y >= 0, z = 0

    cf2 = MOI.ScalarAffineFunction(v, [1.0, -1.0, 0.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf2), MOI.GreaterThan{Float64})
    c2 = MOI.addconstraint!(model, cf2, MOI.GreaterThan(0.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 0

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1, 1, 0] atol=atol rtol=rtol

        @test MOI.canget(model, MOI.ConstraintPrimal(), typeof(c))
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 2 atol=atol rtol=rtol

        if config.duals
            @test MOI.canget(model, MOI.DualStatus(1))
            @test MOI.get(model, MOI.DualStatus(1)) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -1.5 atol=atol rtol=rtol
            @test MOI.get(model, MOI.ConstraintDual(), c2) ≈ 0.5 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc1))
            @test MOI.get(model, MOI.ConstraintDual(), vc1) ≈ 0 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc2))
            @test MOI.get(model, MOI.ConstraintDual(), vc2) ≈ 0 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc3))
            @test MOI.get(model, MOI.ConstraintDual(), vc3) ≈ 1.5 atol=atol rtol=rtol
        end
    end

    if config.query
        @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c2))
        @test MOI.get(model, MOI.ConstraintFunction(), c2) ≈ cf2
    end

    # delete variable x to get :
    # max 2y
    # s.t. y + z == 2
    # - y >= 0
    # y >= 0, z = 0

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 2
    MOI.delete!(model, v[1])
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 1

    if config.query
        @test MOI.get(model, MOI.ConstraintFunction(), c2) ≈ MOI.ScalarAffineFunction([v[2], z], [-1.0, 0.0], 0.0)

        @test MOI.canget(model, MOI.ListOfVariableIndices())
        vrs = MOI.get(model, MOI.ListOfVariableIndices())
        @test vrs == [v[2], z] || vrs == [z, v[2]]
        @test MOI.canget(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()) ≈ MOI.ScalarAffineFunction([v[2], z], [2.0, 0.0], 0.0)
    end
end

# addvariable! (one by one)
function linear2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # Min -x
    # s.t. x + y <= 1
    # x, y >= 0

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}),(MOI.SingleVariable,MOI.GreaterThan{Float64})])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)

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

    cf = MOI.ScalarAffineFunction([x, y], [1.0,1.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf), MOI.LessThan{Float64})
    c = MOI.addconstraint!(model, cf, MOI.LessThan(1.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    vc1 = MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    vc2 = MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 2

    objf = MOI.ScalarAffineFunction([x, y], [-1.0,0.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MinSense

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.ConstraintPrimal(), typeof(c))
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1 atol=atol rtol=rtol

        if config.duals
            @test MOI.canget(model, MOI.DualStatus())
            @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -1 atol=atol rtol=rtol

            # reduced costs
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc1))
            @test MOI.get(model, MOI.ConstraintDual(), vc1) ≈ 0 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(vc2))
            @test MOI.get(model, MOI.ConstraintDual(), vc2) ≈ 1 atol=atol rtol=rtol
        end
    end
end

# Issue #40 from Gurobi.jl
function linear3test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # min  x
    # s.t. x >= 0
    #      x >= 3

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    cf = MOI.ScalarAffineFunction([x], [1.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf), MOI.GreaterThan{Float64})
    MOI.addconstraint!(model, cf, MOI.GreaterThan(3.0))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}}()) == 1

    objf = MOI.ScalarAffineFunction([x], [1.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3 atol=atol rtol=rtol
    end

    # max  x
    # s.t. x <= 0
    #      x <= 3

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.LessThan{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.LessThan(0.0))
    cf = MOI.ScalarAffineFunction([x], [1.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf), MOI.LessThan{Float64})
    MOI.addconstraint!(model, cf, MOI.LessThan(3.0))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.LessThan{Float64}}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    objf = MOI.ScalarAffineFunction([x], [1.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0 atol=atol rtol=rtol
    end
end

# Modify GreaterThan and LessThan sets as bounds
function linear4test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.SingleVariable,MOI.GreaterThan{Float64}),(MOI.SingleVariable,MOI.LessThan{Float64})])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)

    # Min  x - y
    # s.t. 0.0 <= x          (c1)
    #             y <= 0.0   (c2)

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x,y], [1.0, -1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    c1 = MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.LessThan{Float64})
    c2 = MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.LessThan(0.0))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    end

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= 0.0
    @test MOI.canmodifyconstraint(model, c1, MOI.GreaterThan{Float64})
    MOI.modifyconstraint!(model, c1, MOI.GreaterThan(100.0))
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    end

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= -100.0
    @test MOI.canmodifyconstraint(model, c2, MOI.LessThan{Float64})
    MOI.modifyconstraint!(model, c2, MOI.LessThan(-100.0))
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 200.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ -100.0 atol=atol rtol=rtol
    end
end

# Change coeffs, del constr, del var
function linear5test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    #@test MOI.get(model, MOI.SupportsDeleteVariable())
    #####################################
    # Start from simple LP
    # Solve it
    # Copy and solve again
    # Chg coeff, solve, change back solve
    # del constr and solve
    # del var and solve

    #   maximize x + y
    #
    #   s.t. 2 x + 1 y <= 4
    #        1 x + 2 y <= 4
    #        x >= 0, y >= 0
    #
    #   solution: x = 1.3333333, y = 1.3333333, objv = 2.66666666

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)

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

    cf1 = MOI.ScalarAffineFunction([x, y], [2.0,1.0], 0.0)
    cf2 = MOI.ScalarAffineFunction([x, y], [1.0,2.0], 0.0)

    @test MOI.canaddconstraint(model, typeof(cf1), MOI.LessThan{Float64})
    c1 = MOI.addconstraint!(model, cf1, MOI.LessThan(4.0))
    @test MOI.canaddconstraint(model, typeof(cf2), MOI.LessThan{Float64})
    c2 = MOI.addconstraint!(model, cf2, MOI.LessThan(4.0))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 2

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    vc1 = MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    vc2 = MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 2

    objf = MOI.ScalarAffineFunction([x, y], [1.0,1.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 8/3 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), [x, y]) ≈ [4/3, 4/3] atol=atol rtol=rtol
    end

    # copy and solve again
    # missing test

    # change coeff
    #   maximize x + y
    #
    #   s.t. 2 x + 3 y <= 4
    #        1 x + 2 y <= 4
    #        x >= 0, y >= 0
    #
    #   solution: x = 2, y = 0, objv = 2

    @test MOI.canmodifyconstraint(model, c1, MOI.ScalarCoefficientChange{Float64})
    MOI.modifyconstraint!(model, c1, MOI.ScalarCoefficientChange(y, 3.0))
    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), [x, y]) ≈ [2.0, 0.0] atol=atol rtol=rtol
    end

    # delconstrs and solve
    #   maximize x + y
    #
    #   s.t. 1 x + 2 y <= 4
    #        x >= 0, y >= 0
    #
    #   solution: x = 4, y = 0, objv = 4
    @test MOI.candelete(model, c1)
    MOI.delete!(model, c1)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 4 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), [x, y]) ≈ [4.0, 0.0] atol=atol rtol=rtol
    end

    # delvars and solve
    #   maximize y
    #
    #   s.t.  2 y <= 4
    #           y >= 0
    #
    #   solution: y = 2, objv = 2
    @test MOI.candelete(model, x)
    MOI.delete!(model, x)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 2.0 atol=atol rtol=rtol
    end
end

# Modify GreaterThan and LessThan sets as linear constraints
function linear6test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}),(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64})])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)

    # Min  x - y
    # s.t. 0.0 <= x          (c1)
    #             y <= 0.0   (c2)

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x,y], [1.0, -1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})
    c1 = MOI.addconstraint!(model, MOI.ScalarAffineFunction([x],[1.0],0.0), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
    c2 = MOI.addconstraint!(model, MOI.ScalarAffineFunction([y],[1.0],0.0), MOI.LessThan(0.0))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    end

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= 0.0
    @test MOI.canmodifyconstraint(model, c1, MOI.GreaterThan{Float64})
    MOI.modifyconstraint!(model, c1, MOI.GreaterThan(100.0))
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    end

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= -100.0
    @test MOI.canmodifyconstraint(model, c2, MOI.LessThan{Float64})
    MOI.modifyconstraint!(model, c2, MOI.LessThan(-100.0))
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 200.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ -100.0 atol=atol rtol=rtol
    end
end

# Modify constants in Nonnegatives and Nonpositives
function linear7test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.VectorAffineFunction{Float64},MOI.Nonpositives),(MOI.VectorAffineFunction{Float64},MOI.Nonpositives)])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)

    # Min  x - y
    # s.t. 0.0 <= x          (c1)
    #             y <= 0.0   (c2)

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x,y], [1.0, -1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    @test MOI.canaddconstraint(model, MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives)
    c1 = MOI.addconstraint!(model, MOI.VectorAffineFunction([1],[x],[1.0],[0.0]), MOI.Nonnegatives(1))
    @test MOI.canaddconstraint(model, MOI.VectorAffineFunction{Float64}, MOI.Nonpositives)
    c2 = MOI.addconstraint!(model, MOI.VectorAffineFunction([1],[y],[1.0],[0.0]), MOI.Nonpositives(1))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    end

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= 0.0
    @test MOI.canmodifyconstraint(model, c1, MOI.VectorConstantChange{Float64})
    MOI.modifyconstraint!(model, c1, MOI.VectorConstantChange([-100.0]))
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    end

    # Min  x - y
    # s.t. 100.0 <= x
    #               y <= -100.0
    @test MOI.canmodifyconstraint(model, c2, MOI.VectorConstantChange{Float64})
    MOI.modifyconstraint!(model, c2, MOI.VectorConstantChange([100.0]))
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 200.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 100.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ -100.0 atol=atol rtol=rtol
    end
end

# infeasible problem
function linear8atest(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # min x
    # s.t. 2x+y <= -1
    # x,y >= 0
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}),(MOI.SingleVariable,MOI.GreaterThan{Float64})])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
    c = MOI.addconstraint!(model, MOI.ScalarAffineFunction([x,y], [2.0,1.0], 0.0), MOI.LessThan(-1.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    bndx = MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    bndy = MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x], [1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.ResultCount())
        if config.infeas_certificates
            # solver returned an infeasibility ray
            @test MOI.get(model, MOI.ResultCount()) >= 1
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
            @test MOI.get(model, MOI.DualStatus()) == MOI.InfeasibilityCertificate
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            cd = MOI.get(model, MOI.ConstraintDual(), c)
            @test cd < -atol
            # TODO: farkas dual on bounds - see #127
            # xd = MOI.get(model, MOI.ConstraintDual(), bndx)
            # yd = MOI.get(model, MOI.ConstraintDual(), bndy)
            # @test xd > atol
            # @test yd > atol
            # @test yd ≈ -cd atol=atol rtol=rtol
            # @test xd ≈ -2cd atol=atol rtol=rtol
        else
            # solver returned nothing
            @test MOI.get(model, MOI.ResultCount()) == 0
            @test MOI.canget(model, MOI.PrimalStatus(1)) == false
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.InfeasibleNoResult ||
                MOI.get(model, MOI.TerminationStatus()) == MOI.InfeasibleOrUnbounded
        end
    end
end

# unbounded problem
function linear8btest(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # min -x-y
    # s.t. -x+2y <= 0
    # x,y >= 0
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}),(MOI.SingleVariable,MOI.GreaterThan{Float64})])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
    MOI.addconstraint!(model, MOI.ScalarAffineFunction([x,y], [-1.0,2.0], 0.0), MOI.LessThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x, y], [-1.0, -1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.ResultCount())
        if config.infeas_certificates
            # solver returned an unbounded ray
            @test MOI.get(model, MOI.ResultCount()) >= 1
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.InfeasibilityCertificate
        else
            # solver returned nothing
            @test MOI.get(model, MOI.ResultCount()) == 0
            @test MOI.canget(model, MOI.PrimalStatus(1)) == false
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.UnboundedNoResult ||
                MOI.get(model, MOI.TerminationStatus()) == MOI.InfeasibleOrUnbounded
        end
    end
end

# unbounded problem with unique ray
function linear8ctest(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # min -x-y
    # s.t. x-y == 0
    # x,y >= 0
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}),(MOI.SingleVariable,MOI.GreaterThan{Float64})])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})
    MOI.addconstraint!(model, MOI.ScalarAffineFunction([x,y], [1.0,-1.0], 0.0), MOI.EqualTo(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),MOI.ScalarAffineFunction([x, y], [-1.0, -1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.ResultCount())
        if config.infeas_certificates
            # solver returned an unbounded ray
            @test MOI.get(model, MOI.ResultCount()) >= 1
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.InfeasibilityCertificate
            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            ray = MOI.get(model, MOI.VariablePrimal(), [x,y])
            @test ray[1] ≈ ray[2] atol=atol rtol=rtol

        else
            # solver returned nothing
            @test MOI.get(model, MOI.ResultCount()) == 0
            @test MOI.canget(model, MOI.PrimalStatus(1)) == false
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.UnboundedNoResult ||
                MOI.get(model, MOI.TerminationStatus()) == MOI.InfeasibleOrUnbounded
        end
    end
end

# addconstraints
function linear9test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    #   maximize 1000 x + 350 y
    #
    #       s.t.                x >= 30
    #                           y >= 0
    #                 x -   1.5 y >= 0
    #            12   x +   8   y <= 1000
    #            1000 x + 300   y <= 70000
    #
    #   solution: (59.0909, 36.3636)
    #   objv: 71818.1818
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64},
    #    [
    #        (MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}),
    #        (MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}),
    #        (MOI.SingleVariable,MOI.GreaterThan{Float64})
    #    ]
    #)

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)

    MOI.addconstraints!(model,
        [MOI.SingleVariable(x), MOI.SingleVariable(y)],
        [MOI.GreaterThan(30.0), MOI.GreaterThan(0.0)]
    )

    MOI.addconstraints!(model,
        [MOI.ScalarAffineFunction([x, y], [1.0, -1.5], 0.0)],
        [MOI.GreaterThan(0.0)]
    )

    MOI.addconstraints!(model,
        [
            MOI.ScalarAffineFunction([x, y], [12.0, 8.0], 0.0),
            MOI.ScalarAffineFunction([x, y], [1_000.0, 300.0], 0.0)
        ],
        [
            MOI.LessThan(1_000.0),
            MOI.LessThan(70_000.0)
        ]
    )

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
                      MOI.ScalarAffineFunction([x, y], [1_000.0, 350.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 79e4/11 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 650/11 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 400/11 atol=atol rtol=rtol
    end
end

# ranged constraints
function linear10test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    #   maximize x + y
    #
    #       s.t.  5 <= x + y <= 10
    #                  x,  y >= 0
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64},
    #    [
    #        (MOI.ScalarAffineFunction{Float64},MOI.Interval{Float64}),
    #        (MOI.SingleVariable,MOI.GreaterThan{Float64})
    #    ]
    #)

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)

    MOI.addconstraints!(model,
        [MOI.SingleVariable(x), MOI.SingleVariable(y)],
        [MOI.GreaterThan(0.0), MOI.GreaterThan(0.0)]
    )

    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
    c = MOI.addconstraint!(model, MOI.ScalarAffineFunction([x,y], [1.0, 1.0], 0.0), MOI.Interval(5.0, 10.0))

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x, y], [1.0, 1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 10.0 atol=atol rtol=rtol

        if config.duals
            @test MOI.get(model, MOI.ResultCount()) >= 1
            @test MOI.canget(model, MOI.DualStatus())
            @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -1 atol=atol rtol=rtol
        end
    end

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x, y], [1.0, 1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 5.0 atol=atol rtol=rtol

        if config.duals
            @test MOI.get(model, MOI.ResultCount()) >= 1
            @test MOI.canget(model, MOI.DualStatus())
            @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            @test MOI.get(model, MOI.ConstraintDual(), c) ≈ 1 atol=atol rtol=rtol
        end
    end

    @test MOI.canmodifyconstraint(model, c, MOI.Interval{Float64})
    MOI.modifyconstraint!(model, c, MOI.Interval(2.0, 12.0))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
    end

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x, y], [1.0, 1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 12.0 atol=atol rtol=rtol
    end
end

# changing constraint sense
function linear11test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # simple 2 variable, 1 constraint problem
    # min x + y
    # st   x + y >= 1
    #      x + y >= 2
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64},
    #    [
    #        (MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}),
    #        (MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64})
    #    ]
    #)

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    v = MOI.addvariables!(model, 2)

    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})
    c1 = MOI.addconstraint!(model, MOI.ScalarAffineFunction(v, [1.0,1.0], 0.0), MOI.GreaterThan(1.0))
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})
    c2 = MOI.addconstraint!(model, MOI.ScalarAffineFunction(v, [1.0,1.0], 0.0), MOI.GreaterThan(2.0))

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(v, [1.0,1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
    end

    @test MOI.cantransformconstraint(model, c2, MOI.LessThan{Float64})
    c3 = MOI.transformconstraint!(model, c2, MOI.LessThan(2.0))

    @test isa(c3, MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}})
    @test MOI.isvalid(model, c2) == false
    @test MOI.isvalid(model, c3) == true

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol
    end
end

# infeasible problem with 2 linear constraints
function linear12test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # min x
    # s.t. 2x-3y <= -7
    #      y <= 2
    # x,y >= 0
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}),(MOI.SingleVariable,MOI.GreaterThan{Float64})])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    x = MOI.addvariable!(model)
    MOI.canaddvariable(model)
    y = MOI.addvariable!(model)
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
    c1 = MOI.addconstraint!(model, MOI.ScalarAffineFunction([x,y], [2.0,-3.0], 0.0), MOI.LessThan(-7.0))
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
    c2 = MOI.addconstraint!(model, MOI.ScalarAffineFunction([y], [1.0], 0.0), MOI.LessThan(2.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    bndx = MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    bndy = MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x], [1.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.ResultCount())
        if config.infeas_certificates
            # solver returned an infeasibility ray
            @test MOI.get(model, MOI.ResultCount()) >= 1
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
            @test MOI.get(model, MOI.DualStatus()) == MOI.InfeasibilityCertificate
            @test MOI.canget(model, MOI.ConstraintDual(), typeof(c1))
            cd1 = MOI.get(model, MOI.ConstraintDual(), c1)
            cd2 = MOI.get(model, MOI.ConstraintDual(), c2)
            bndxd = MOI.get(model, MOI.ConstraintDual(), bndx)
            bndyd = MOI.get(model, MOI.ConstraintDual(), bndy)
            @test cd1 < - atol
            @test cd2 < - atol
            @test - 3 * cd1 + cd2 ≈ -bndyd atol=atol rtol=rtol
            @test 2 * cd1 ≈ -bndxd atol=atol rtol=rtol
            @test -7 * cd1 + 2 * cd2 > atol
        else
            # solver returned nothing
            @test MOI.get(model, MOI.ResultCount()) == 0
            @test MOI.canget(model, MOI.PrimalStatus(1)) == false
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.InfeasibleNoResult ||
                MOI.get(model, MOI.TerminationStatus()) == MOI.InfeasibleOrUnbounded
        end
    end
end

const contlineartests = Dict("linear1" => linear1test,
                              "linear2" => linear2test,
                              "linear3" => linear3test,
                              "linear4" => linear4test,
                              "linear5" => linear5test,
                              "linear6" => linear6test,
                              "linear7" => linear7test,
                              "linear8a" => linear8atest,
                              "linear8b" => linear8btest,
                              "linear8c" => linear8ctest,
                              "linear9" => linear9test,
                              "linear10" => linear10test,
                              "linear11" => linear11test,
                              "linear12" => linear12test)

@moitestset contlinear


text:
### Vector Policy ###
# maintained by @zsunberg and @etotheipluspi

"""
A generic MDP policy that consists of a vector of actions. The entry at `state_index(mdp, s)` is the action that will be taken in state `s`.
"""
mutable struct VectorPolicy{S,A} <: Policy
    mdp::MDP{S,A}
    act::Vector{A}
end

action(p::VectorPolicy, s) = p.act[state_index(p.mdp, s)]
action(p::VectorPolicy, s, a) = action(p, s)

"""
Solver for VectorPolicy. Doesn't do any computation - just sets the action vector.
"""
mutable struct VectorSolver{A}
    act::Vector{A}
end

create_policy(s::VectorSolver{A}, mdp::MDP{S,A}) where {S,A} = VectorPolicy(mdp, Array{A}(0))

function solve(s::VectorSolver, mdp::MDP, p::VectorPolicy=create_policy(s,mdp))
    p.mdp = mdp
    p.act = s.act
    return p
end

"""
A generic MDP policy that consists of a value table. The entry at `state_index(mdp, s)` is the action that will be taken in state `s`.
"""
mutable struct ValuePolicy{A, P <: Union{MDP, POMDP}} <: Policy
    mdp::P
    value_table::Matrix{Float64}
    act::Vector{A}
end
function ValuePolicy(mdp::P) where {P <: Union{MDP, POMDP}}
    A = action_type(mdp)
    return ValuePolicy{A, P}(mdp, zeros(n_states(mdp), n_actions(mdp)), ordered_actions(mdp))
end

action(p::ValuePolicy, s) = p.act[indmax(p.value_table[state_index(p.mdp, s),:])]

# @POMDP_require ValuePolicy(mdp::Union{MDP, POMDP}) begin
#     M = typeof(mdp)
#     @req n_states(::M)
#     @req n_actions(::M)
#     @subreq ordered_actions(mdp)
# end
#
# @POMDP_require action(p::ValuePolicy, s) begin
#     @req state_index(::typeof(p.mdp), ::typeof(s))
# end


text:
abstract type IntegratorCache{DT,D} end

abstract type ODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type DAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end

IntegratorCache(params::Parameters) = error("IntegratorCache(params) not implemented for ", typeof(params))

create_integrator_cache(integrator::Integrator) = error("create_integrator_cache()! not implemented for ", typeof(integrator))

initialize!(::Integrator, ::AtomicSolution) = nothing

integrate_step!(integrator::Integrator, ::AtomicSolution) = error("integrate_step()! not implemented for ", typeof(integrator))


text:
using Test
using Random
using TiledViews

@testset "testing array access" begin
    a = TiledView(reshape(1:49,(7,7)), (4, 4),(0, 0));
    @test size(a) == (4,4,2,2)
    @test a[1,1,1,1] == 0
    @test a[4,4,1,1] == 17
end

@testset "Check errors" begin
    a = TiledView(reshape(1:49,(7,7)), (4, 4),(0, 0));
    @test_throws BoundsError a[0,1,1,1]
    @test_throws BoundsError a[5,1,1,1]
    @test_throws BoundsError a[4,0,1,1]
    @test_throws BoundsError a[4,5,1,1]
    @test_throws BoundsError a[4,1,0,1]
    @test_throws BoundsError a[4,1,4,1]
    @test_throws BoundsError a[4,1,1,0]
    @test_throws BoundsError a[4,1,1,4]
end

@testset "Test own size method" begin
    q = ones(7,7) .+ 0.0
    a = TiledView(q, (5, 4),(3, 2));
    a[:,:,:,:] .= 77;
    @test !any(q .!= 77.0)
end

@testset "Test tile_centers" begin
    q = rand(100,101);
    a=TiledView(q, (80,91), (0,0));
    @test [(-20, -10) (-20, 81);  (60, -10) (60, 81)] == tile_centers(a)
end

#=
@testset "similar" begin
    function compare_sizes(a)
        b=similar(a);
        #print(size(a));print(size(b));size(b.parent)
        @test all(size(a) .== size(b))
        @test all(size(a.parent) .<= size(b.parent))
        c=similar(a,Int64,size(b).+1)
        @test all(size(c) .== size(a) .+1)
        @test all(size(c.parent) .> size(a.parent))
    end
    for n in 1:40
        dims = rand(1:5);
        sz = Tuple(rand(3:7,dims));
        ts = Tuple(rand(1:4,dims));
        overlap = Tuple(rand(0:2,dims));
        overlap = Tuple(((overlap[n] >= ts[n]) ? ts[n]-1 : overlap[n] for n in 1:length(ts)));
        keep_ctr = rand(Bool);
        a = TiledView(ones(sz...), ts, overlap, keep_center=keep_ctr);
        compare_sizes(a);
    end
end
=#
@testset "zeros_like" begin
    q = rand(100,101);
    a=TiledView(q, (80,91), (0,0), keep_center=false);
    w = TiledViews.zeros_like(a)
    @test size(w) == size(a)
end

@testset "similar" begin
    q = rand(100,101);
    a=TiledView(q, (80,91), (0,0));
    @test size(similar(a)) == size(a)
    @test eltype(similar(a,Int32)) == Int32
    @test size(similar(a,Int32,(30,40,1))) == (30,40,1)
end

@testset "TiledWindowView" begin
    dest = ones(512,512)
    tile_size = (65,65)
    tiles, win = TiledWindowView(dest, tile_size, rel_overlap=(1.0,1.0));
    to_write = 0.2 .* win .* collect(tiles[:,:,:,:]);
    #tiles[:,:,:,:] .+= to_write[:,:,:,:] ;  # why can one NOT use [:,:,:,:]?
    tiles .+= to_write ;  # why can one NOT use [:,:,:,:]?
    @test sum(abs.(dest[50:end-50,50:end-50] .- 1.2)) < 1e-10
    win2, mynorm = get_window(tiles, get_norm=true, verbose=true);
    @test all(win .== win2)
    tiles, win = TiledWindowView(dest, tile_size, rel_overlap=(1.0,1.0), get_norm=true);
    @test all(win .== win2)
end

@testset "sub-indexing" begin
    q = rand(100,101);
    a=TiledView(q, (80,91), (0,0));
    @test size(a[:,3:14,1:2,1]) == (80,12,2)
end

@testset "eachtile" begin
    q = zeros(100,101);
    a=TiledView(q, (80,91), (0,0));
    for (t,tn,tc,tc2) in zip(eachtile(a),eachtilenumber(a),eachtilerelpos(a),eachtilerelpos(a,2.0))
        t[:] .+= tn[1];
    end
    @test q[end,end] == 2
end

@testset "tiled_processing" begin
    q = ones(10,10,10);
    a = TiledView(q, (4,4,4), (2,2,2));
    fct(a)=sin.(a)
    res = tiled_processing(q, fct, (4,4,4), (2,2,2))
    @test maximum(res.parent) ≈ sin(1)
end

return


text:
using GLMakie

points = [Point2f0(cos(t), sin(t)) for t in LinRange(0, 2pi, 20)]
colors = 1:20
figure, axis, scatterobject = scatter(points, color = colors, markersize = 15)
figure

max_stars_repo_name: JMSoukup/github-demo

max_issues_repo_name: JMSoukup/github-demo

max_forks_repo_name: JMSoukup/github-demo


text:
"""
    evalcompositelogisticprobit(x::T, a, b)::T where T <: Real

evaluates the map from (0,1) to (0,1).
returns 1/(1 + exp(-a*(log(x/(1-x))-b)))

    evalcompositelogisticprobit(x::T, a, b, lb, ub)::T where T <: Real
evaluates the map from (lb, ub) to (lb, ub), by mapping (lb,ub) to (0,1), run evalcompositelogisticprobit(), then map (0,1) to (lb,ub).
"""
function evalcompositelogisticprobit(x::T, a, b)::T where T <: Real

    return 1/(1 + exp(-a*(log(x/(1-x))-b)))
end

function evalcompositelogisticprobit(x_inp::T, a, b, lb, ub)::T where T <: Real

    x = convertcompactdomain(x_inp, lb, ub, zero(T), one(T))
    y = evalcompositelogisticprobit(x, a, b)
    y_out = convertcompactdomain(y, zero(T), one(T), lb, ub)

    return y_out
end

function evalinversecompositelogisticprobit(y::T, a, b)::T where T <: Real
    return exp(b)/(exp(b) + (-1 + 1/y)^(1/a))
end

function evalinversecompositelogisticprobit(y_inp::T, a, b, lb, ub)::T where T <: Real

    y = convertcompactdomain(y_inp, lb, ub, zero(T), one(T))
    x = evalinversecompositelogisticprobit(y, a, b)
    x_out = convertcompactdomain(x, zero(T), one(T), lb, ub)

    return x_out
end

function eval1Dnumericalinverse(f::Function,
    y::T,
    x0::T,
    a::T,
    b::T,
    max_iters::Int) where T <: Real

    @assert a < b

    obj_func = xx->((f(xx[1])-y)^2)::T

    op = Optim.Options( iterations = max_iters,
                         store_trace = false,
                         show_trace = false)

    results = Optim.optimize(   obj_func,
                                [x0],
                                Optim.NewtonTrustRegion(),
                                op)

    x_star = results.minimizer
    x_out = clamp(x_star[1], a, b)

    return x_out, results
end


text:
module DataTypes

using StaticArrays
using JuLIP: vecs, mat, JVec, energy, forces, virial
import IPFitting: vec_obs, devec_obs, eval_obs, weighthook, Dat, err_weighthook

export ENERGY, FORCES, VIRIAL

# some convenience functions to dispatch string arguments
vec_obs(s::AbstractString, args...) = vec_obs(Val(Symbol(s)), args...)
devec_obs(s::AbstractString, args...) = devec_obs(Val(Symbol(s)), args...)
eval_obs(s::AbstractString, args...) = eval_obs(Val(Symbol(s)), args...)
weighthook(s::AbstractString, args...) = weighthook(Val(Symbol(s)), args...)

# ------------------- DEFAULTS ------------------

vec_obs(::Val, x::AbstractVector{<: AbstractFloat}) = collect(x)
devec_obs(::Val,  x::AbstractVector) = x

# ------------------- ENERGY ------------------

const ENERGY = "E"
const ValE = Val{:E}
vec_obs(::ValE, E::Real) = [E]
vec_obs(::ValE, E::Vector{<: Real}) = E
devec_obs(::ValE, x::AbstractVector) = ((@assert length(x) == 1); x[1])
eval_obs(::ValE, B, dat::Dat) = energy(B, dat.at)
weighthook(::ValE, d::Dat) = 1.0 / sqrt(length(d.at))
err_weighthook(::ValE, d::Dat) = 1.0 / length(d.at)

# ------------------- FORCES ------------------

const FORCES = "F"
const ValF = Val{:F}
vec_obs(v::ValF, F::AbstractVector{<:JVec}) = vec_obs(v, mat(F))
vec_obs(::ValF, F::AbstractMatrix) = vec(F)
devec_obs(::ValF, x::AbstractVector) = vecs(reshape(x, 3, :))
eval_obs(::ValF, B, dat::Dat) = forces(B, dat.at)

function vec_obs(valF::ValF, F::Vector{<: Vector})
   nbasis = length(F)
   nat = length(F[1])
   Fmat = zeros(3*nat, nbasis)
   for ib = 1:nbasis
      Fmat[:, ib] .= vec_obs(valF, F[ib])
   end
   return Fmat
end

# ------------------- VIRIAL ------------------
# using Voigt convention for vectorising  symmetric 3 x 3 matrix
#  1  6  5
#  6  2  4
#  5  4  3
const _IV = [1,5,9,6,3,2]
const _IVst = SVector(1,5,9,6,3,2)
const VIRIAL = "V"
const ValV = Val{:V}
vec_obs(::ValV, v::AbstractVector{<: Real}) = (@assert length(v) == 6; collect(v))
vec_obs(::ValV, V::AbstractMatrix) = (@assert size(V) == (3,3); V[_IV])
devec_obs(::ValV, x::AbstractVector) =
   SMatrix{3,3}(x[1], x[6], x[5], x[6], x[2], x[4], x[5], x[4], x[3])
eval_obs(::ValV, B, dat::Dat) = virial(B, dat.at)
weighthook(::ValV, d::Dat) = 1.0 / sqrt(length(d.at))
err_weighthook(::ValV, d::Dat) = 1.0 / length(d.at)

function vec_obs(valV::ValV, V::AbstractVector{<: AbstractMatrix})
   nbasis = length(V)
   Vmat = zeros(6, nbasis)
   for ib = 1:nbasis
      Vmat[:, ib] .= vec_obs(valV, V[ib])
   end
   return Vmat
end

end


text:
#include("coordinates.jl")
#using .Coordinatesm

@enum ResourceName res_type_a=1 res_type_b
res_type_energy=[20,50]

mutable struct ResourceType
    name::ResourceName
    energy::Int64
    function ResourceType(name::ResourceName)
        new(name,res_type_energy[Int(name)])
    end
end

mutable struct Resource
    id::Int64
    type::ResourceType
    position::Coordinates
    free_energy::Int64
end

function create_random_resource(id::Int,possible_coordinates::PossibleCoordinates)
    res_type=ResourceType(rand(instances(ResourceName)))
    Resource(id,res_type,get_random_coordinates(possible_coordinates),res_type.energy)
end


text:
# core definition of the AbstractGPUArray type

export AbstractGPUArray

"""
    AbstractGPUArray{T, N} <: DenseArray{T, N}

Supertype for `N`-dimensional GPU arrays (or array-like types) with elements of type `T`.
Instances of this type are expected to live on the host, see [`AbstractDeviceArray`](@ref)
for device-side objects.
"""
abstract type AbstractGPUArray{T, N} <: DenseArray{T, N} end

const AbstractGPUVector{T} = AbstractGPUArray{T, 1}
const AbstractGPUMatrix{T} = AbstractGPUArray{T, 2}
const AbstractGPUVecOrMat{T} = Union{AbstractGPUArray{T, 1}, AbstractGPUArray{T, 2}}

# convenience aliases for working with wrapped arrays

const WrappedGPUArray{T,N} = WrappedArray{T,N,AbstractGPUArray,AbstractGPUArray{T,N}}

const AnyGPUArray{T,N} = Union{AbstractGPUArray{T,N}, WrappedGPUArray{T,N}}

# input/output

## serialization

using Serialization: AbstractSerializer, serialize_type

function Serialization.serialize(s::AbstractSerializer, t::T) where T <: AbstractGPUArray
    serialize_type(s, T)
    serialize(s, Array(t))
end

function Serialization.deserialize(s::AbstractSerializer, ::Type{T}) where T <: AbstractGPUArray
    A = deserialize(s)
    T(A)
end

## showing

struct ToArray end
Adapt.adapt_storage(::ToArray, xs::AbstractGPUArray) = convert(Array, xs)

# display
Base.print_array(io::IO, X::AnyGPUArray) =
    Base.print_array(io, adapt(ToArray(), X))

# show
Base._show_nonempty(io::IO, X::AnyGPUArray, prefix::String) =
    Base._show_nonempty(io, adapt(ToArray(), X), prefix)
Base._show_empty(io::IO, X::AnyGPUArray) =
    Base._show_empty(io, adapt(ToArray(), X))
Base.show_vector(io::IO, v::AnyGPUArray, args...) =
    Base.show_vector(io, adapt(ToArray(), v), args...)

## collect to CPU (discarding wrapper type)

collect_to_cpu(xs::AbstractArray) = collect(adapt(ToArray(), xs))
Base.collect(X::AnyGPUArray) = collect_to_cpu(X)

# memory copying

## basic linear copies of identically-typed memory

# expects the GPU array type to have linear `copyto!` methods (i.e. accepting an integer
# offset and length) from and to CPU arrays and between GPU arrays.

for (D, S) in ((AnyGPUArray, Array),
               (Array, AnyGPUArray),
               (AnyGPUArray, AnyGPUArray))
    @eval begin
        function Base.copyto!(dest::$D{<:Any, N}, rdest::UnitRange,
                              src::$S{<:Any, N}, ssrc::UnitRange) where {N}
            drange = CartesianIndices((rdest,))
            srange = CartesianIndices((ssrc,))
            copyto!(dest, drange, src, srange)
        end

        function Base.copyto!(dest::$D, d_range::CartesianIndices{1},
                              src::$S, s_range::CartesianIndices{1})
            len = length(d_range)
            if length(s_range) != len
                throw(ArgumentError("Copy range needs same length. Found: dest: $len, src: $(length(s_range))"))
            end
            len == 0 && return dest
            d_offset = first(d_range)[1]
            s_offset = first(s_range)[1]
            copyto!(dest, d_offset, src, s_offset, len)
        end

        Base.copyto!(dest::$D, src::$S) = copyto!(dest, 1, src, 1, length(src))
    end
end

# kernel-based variant for copying between wrapped GPU arrays

function linear_copy_kernel!(ctx::AbstractKernelContext, dest, dstart, src, sstart, n)
    i = linear_index(ctx)-1
    if i < n
        @inbounds dest[dstart+i] = src[sstart+i]
    end
    return
end

function Base.copyto!(dest::AnyGPUArray, dstart::Integer,
                      src::AnyGPUArray, sstart::Integer, n::Integer)
    n == 0 && return dest
    n < 0 && throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative")))
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    (checkbounds(Bool, destinds, dstart) && checkbounds(Bool, destinds, dstart+n-1)) || throw(BoundsError(dest, dstart:dstart+n-1))
    (checkbounds(Bool, srcinds, sstart)  && checkbounds(Bool, srcinds, sstart+n-1))  || throw(BoundsError(src,  sstart:sstart+n-1))

    gpu_call(linear_copy_kernel!,
             dest, dstart, src, sstart, n;
             elements=n)
    return dest
end

# variants that materialize the GPU wrapper before copying from or to the CPU

function Base.copyto!(dest::Array, dstart::Integer,
                      src::WrappedGPUArray, sstart::Integer, n::Integer)
    n == 0 && return dest
    temp = similar(src, n)
    copyto!(temp, 1, src, sstart, n)
    copyto!(dest, dstart, temp, 1, n)
    return dest
end

function Base.copyto!(dest::WrappedGPUArray, dstart::Integer,
                      src::Array, sstart::Integer, n::Integer)
    n == 0 && return dest
    temp = similar(dest, n)
    copyto!(temp, 1, src, sstart, n)
    copyto!(dest, dstart, temp, 1, n)
    return dest
end

# variants that converts values on the CPU when there's a type mismatch
#
# we prefer to convert on the CPU where there's typically more memory / less memory pressure
# to quickly perform these very lightweight conversions

function Base.copyto!(dest::Array{T}, dstart::Integer,
                      src::AnyGPUArray{U}, sstart::Integer,
                      n::Integer) where {T,U}
    n == 0 && return dest
    temp = Vector{U}(undef, n)
    copyto!(temp, 1, src, sstart, n)
    copyto!(dest, dstart, temp, 1, n)
    return dest
end

function Base.copyto!(dest::AnyGPUArray{T}, dstart::Integer,
                      src::Array{U}, sstart::Integer, n::Integer) where {T,U}
    n == 0 && return dest
    temp = Vector{T}(undef, n)
    copyto!(temp, 1, src, sstart, n)
    copyto!(dest, dstart, temp, 1, n)
    return dest
end

## generalized blocks of heterogeneous memory

function cartesian_copy_kernel!(ctx::AbstractKernelContext, dest, dest_offsets, src, src_offsets, shape, length)
    i = linear_index(ctx)
    if i <= length
        idx = CartesianIndices(shape)[i]
        @inbounds dest[idx + dest_offsets] = src[idx + src_offsets]
    end
    return
end

function Base.copyto!(dest::AnyGPUArray{<:Any, N}, destcrange::CartesianIndices{N},
                      src::AnyGPUArray{<:Any, N}, srccrange::CartesianIndices{N}) where {N}
    shape = size(destcrange)
    if shape != size(srccrange)
        throw(ArgumentError("Ranges don't match their size. Found: $shape, $(size(srccrange))"))
    end
    len = length(destcrange)
    len == 0 && return dest

    dest_offsets = first(destcrange) - oneunit(CartesianIndex{N})
    src_offsets = first(srccrange) - oneunit(CartesianIndex{N})
    gpu_call(cartesian_copy_kernel!,
             dest, dest_offsets, src, src_offsets, shape, len;
             elements=len)
    dest
end

for (dstTyp, srcTyp) in (AbstractGPUArray=>Array, Array=>AbstractGPUArray)
    @eval function Base.copyto!(dst::$dstTyp{T,N}, dstrange::CartesianIndices{N},
                                src::$srcTyp{T,N}, srcrange::CartesianIndices{N}) where {T,N}
        isempty(dstrange) && return dst
        if size(dstrange) != size(srcrange)
            throw(ArgumentError("source and destination must have same size (got $(size(srcrange)) and $(size(dstrange)))"))
        end

        # figure out how many dimensions of the Cartesian ranges map onto contiguous memory
        # in both source and destination. we will copy these one by one as linear ranges.
        contiguous_dims = 1
        for dim in 2:N
            # a slice is broken up if the previous dimension didn't cover the entire range
            if axes(src, dim-1) == axes(srcrange, dim-1) &&
            axes(dst, dim-1) == axes(dstrange, dim-1)
                contiguous_dims = dim
            else
                break
            end
        end

        m = prod(size(dstrange)[1:contiguous_dims])       # inner, contiguous length
        n = prod(size(dstrange)[contiguous_dims+1:end])   # outer non-contiguous length
        @assert m*n == length(srcrange) == length(dstrange)

        # copy linear slices
        for i in 1:m:m*n
            srcoff = LinearIndices(src)[srcrange[i]]
            dstoff = LinearIndices(dst)[dstrange[i]]
            # TODO: Use asynchronous memory copies
            copyto!(dst, dstoff, src, srcoff, m)
        end

        dst
    end
end

## other

Base.copy(x::AbstractGPUArray) = error("Not implemented") # COV_EXCL_LINE

Base.deepcopy(x::AbstractGPUArray) = copy(x)

# filtering

# TODO: filter!

# revert of JuliaLang/julia#31929
Base.filter(f, As::AbstractGPUArray) = As[map(f, As)::AbstractGPUArray{Bool}]


text:
@test norm(integrate(x -> 1 / x, 1.0, 2.0) - log(2)) < 10e-8
@test norm(integrate(x -> -sin(x), 0.0, pi) - (cos(pi) - cos(0.0))) < 10e-8

r = integrate(x -> 1, 0, 1)
@test norm(r - 1) < 10e-8

r = integrate(x -> x, 0, 1)
@test norm(r - 0.5) < 10e-8

r = integrate(x -> x * x, 0, 1)
@test norm(r - 1 / 3) < 10e-8

r = integrate(sin, 0, pi)
@test norm(r - 2) < 10e-8

r = integrate(cos, 0, pi)
@test norm(r - 0) < 10e-8

r = integrate(x -> sin(x)^2 + sin(x)^2, 0, pi)
@test norm(r - pi) < 10e-8

# Nice example, but requires Distributions
# require("Distributions")
# using Distributions
# r = integrate(x -> pdf(Normal(0.0, 1.0), x), -10, 10)
# @test norm(1 - r) < 10e-8

r = integrate(x -> 1 / x, 0.01, 1)
@test norm(r - 4.60517) < 10e-7

r = integrate(x -> cos(x)^8, 0, 2 * pi)
@test norm(r - 35 * pi / 64) < 10e-7

r = integrate(x -> sin(x) - sin(x^2) + sin(x^3), pi, 2 * pi)
@test norm(r - (-1.830467)) < 10e-7

# Monte Carlo integration tests
r = integrate(x -> sin(x), 0, pi, :monte_carlo)
@test norm(r - 2) < 10e-1


text:
using AutomotiveDrivingModels
using NGSIM
using HDF5

include("RootDir.jl")
const TRAIN_TEST_SPLIT_FILEPATH = joinpath(ROOT_FILEPATH, "data", "NGSIM_train_test_split.h5")

export
    load_assignment,
    load_evaldata

function load_assignment(filepath::AbstractString = TRAIN_TEST_SPLIT_FILEPATH)
    h5read(filepath, "data/assignment")
end

function load_evaldata(filepath::AbstractString = TRAIN_TEST_SPLIT_FILEPATH)
    trajdata_indeces = h5read(filepath, "data/trajdata_indeces")
    egoids = h5read(filepath, "data/egoids")
    frame_lo = h5read(filepath, "data/frame_lo")
    frame_hi = h5read(filepath, "data/frame_hi")

    trajdatas = map(i->load_trajdata(i), 1:length(TRAJDATA_PATHS))
    segments = map(i->TrajdataSegment(trajdata_indeces[i], egoids[i], frame_lo[i], frame_hi[i]), 1:length(trajdata_indeces))

    EvaluationData(trajdatas, segments)
end

max_stars_repo_name: yoshihikokuwahara/gail-driver

max_issues_repo_name: yoshihikokuwahara/gail-driver

max_forks_repo_name: yoshihikokuwahara/gail-driver


text:
module Benchmark
using StaticArrays, PyCall
using AbstractTrees
using Parameters, Combinatorics

import ..GenericPara

const DI, EX = 1, 2
const INL, OUTL, INR, OUTR = 1, 2, 3, 4
# orginal diagrams T, U, S; particle-hole counterterm Ts, Us; and their counterterm Tc, Uc, Sc, Tsc, Usc
const I, T, U, S, Ts, Us, Ic, Tc, Uc, Sc, Tsc, Usc = 1:12
# const ChanName = ["I", "T", "U", "S", "Ts", "Us", "Ic", "Tc", "Uc", "Sc", "Tsc", "Usc"]
const SymFactor = [1.0, -1.0, 1.0, -0.5, +1.0, -1.0]

# const Fchan = [I, U, S, Ts, Us, Ic, Uc, Sc, Tsc, Usc]
# const Vchan = [I, T, U, Ts, Us, Ic, Tc, Uc, Sc, Tsc, Usc]
# const Allchan = [I, T, U, S, Ts, Us, Ic, Tc, Uc, Sc, Tsc, Usc]

include("vertex4.jl")
include("vertex4_eval.jl")
include("vertex4_io.jl")
include("diagram_count.jl")
end


text:
module BitcodeRunner

export bitcode_library

using LLVM

# Internal convertions used for arguments to `llvmcall`.
myconvert(T, x) = Base.cconvert(T, x)
myconvert(T::Type{Base.Ptr{V}}, x) where {V} = Base.unsafe_convert(T, x)

"""
    bitcode_library(bitcodefile, mod::Module = Module())

Reads in LLVM bitcode from the file named `bitcodefile`. Populates `mod` with Julia
functions that map to every function defined in the bitcode file. This works by defining
Julia functions in `mod` that use `Base.llvmcall` to run the LLVM code. `mod` optionally
defaults to an anonymous Module.

Returns the Module with the loaded functions.
"""
function bitcode_library(bitcodefile, mod::Module = Module())
    ctx = LLVM.Context(convert(LLVM.API.LLVMContextRef, cglobal(:jl_LLVMContext, Void)))
    llvmmod = parse(LLVM.Module, read(bitcodefile), ctx)
    fns = functions(llvmmod)
    exs = Any[]
    for f in fns
            # @show attr = collect(function_attributes(f))
            # EnumAttribute("internal") in attr && continue
            # EnumAttribute("hidden") in attr && continue
            funname = Symbol(name(f))
        try
            rettypeX, argtypesX = params(f)
            rettype = maptype(rettypeX, mod)
            argtypes = maptype.(argtypesX, mod)
            n = length(argtypes)
            argnames = [Symbol(:x, i) for i in 1:n]
            # push!(exs, :(export $funname))
            # push!(exs, :($funname($(argnames...)) = Base.llvmcall($(LLVM.ref(f)), $rettype, Tuple{$(argtypes...)}, $([:($myconvert($(argtypes[i]), $(argnames[i]))) for i in 1:n]...))))
            # println(:($funname($(argnames...)) = Base.llvmcall($(LLVM.ref(f)), $rettype, Tuple{$(argtypes...)}, $([:($myconvert($(argtypes[i]), $(argnames[i]))) for i in 1:n]...))))
            eval(mod, :(export $funname; $funname($(argnames...)) = Base.llvmcall($(LLVM.ref(f)), $rettype, Tuple{$(argtypes...)}, $([:($myconvert($(argtypes[i]), $(argnames[i]))) for i in 1:n]...))))
            println("$funname")
        catch
            println("***FAILED: $funname")
        end
    end
    # open("t.out", "w") do io; print(io, Expr(:block, exs...)) end
    # eval(mod, Expr(:block, exs...))
    return mod
end

#
# Find the return type and the parameters of the LLVM function `f`.
# `returntype` is an LLVM type.
# `paramtypes` is a Vector of LLVM types, one for each argument.
#
function params(f)
    funtype = LLVM.API.LLVMGetElementType(LLVM.API.LLVMTypeOf(LLVM.ref(f)))
    returntype = LLVMType(LLVM.API.LLVMGetReturnType(funtype))
    nparams = LLVM.API.LLVMCountParamTypes(funtype)
    paramtypesraw = Vector{LLVM.API.LLVMTypeRef}(nparams)
    LLVM.API.LLVMGetParamTypes(funtype, paramtypesraw)
    paramtypes = LLVMType.(paramtypesraw)
    (returntype, paramtypes)
end

#
# Map the LLVM type (i32*) to its matching Julia type (Ptr{Int32}).
#
function maptype(x, m::Module)
    if isa(x, LLVM.VoidType)
        return Void
    end
    if isa(x, LLVM.IntegerType)
        w = width(x)
        return w == 1 ? Bool : eval(Symbol(:Int, w))
    end
    if isa(x, LLVM.FloatingPointType)
        return isa(x, LLVM.LLVMHalf)   ? Float16 :
               isa(x, LLVM.LLVMFloat)  ? Float32 :
               isa(x, LLVM.LLVMDouble) ? Float64 : nothing
    end
    if isa(x, LLVM.PointerType)
        return Ptr{maptype(eltype(x), m)}
    end
    if isa(x, LLVM.StructType)
        typename = Symbol(split(name(x), ".")[2])
        if !isdefined(m, typename)  # Create the type in Module `m` if not already there
            element_types = map(x -> maptype(x, m), elements(x))
            type_definition = :(struct $typename  end)
            type_definition.args[3].args = Any[:($(Symbol(:x,i)) :: $(element_types[i])) for i in 1:length(element_types)]
            eval(m, type_definition)
        end
        return getfield(m, typename)
    end
    if isa(x, LLVM.ArrayType)
        return NTuple{Int32(length(x)), maptype(eltype(x), m)}
    end
    error("Unknown or unsupported LLVM type: $x")
end

end # module