# Connecting a simple first-order solver to solve standard form linear program to JuMP

**URL:** <https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [March 7, 2023, 8:56pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694 "2023-03-07T20:56:35Z")\
**Posts on this page:** 18\
**Page:** 2

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 14, 2023, 11:13pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/23 "2023-08-14T23:13:04Z")

</div>

If you want to model `ScalarAffineFunction-in-GreaterThan` with `VariableIndex`-in-`Interval` constraints, then you’ll need something like the following. (It isn’t sufficient to just define `supports_constraint`, you also need to setup the data structure accordingly.)

```julia
import MathOptInterface as MOI

MOI.Utilities.@product_of_sets(RHS, MOI.GreaterThan{T})

const OptimizerCache = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            Int,
            MOI.Utilities.OneBasedIndexing,
        },
        Vector{Float64},
        RHS{Float64},
    },
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    x_primal::Dict{MOI.VariableIndex,Float64}
    termination_status::MOI.TerminationStatusCode

    function Optimizer()
        return new(Dict{MOI.VariableIndex,Float64}(), MOI.OPTIMIZE_NOT_CALLED)
    end
end

function MOI.is_empty(model::Optimizer)
    return isempty(model.x_primal) &&
        model.termination_status == MOI.OPTIMIZE_NOT_CALLED
end

function MOI.empty!(model::Optimizer)
    empty!(model.x_primal)
    model.termination_status = MOI.OPTIMIZE_NOT_CALLED
    return
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.GreaterThan{Float64}},
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.EqualTo{Float64},
            MOI.Interval{Float64},
        },
    }
)
    return true
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = OptimizerCache()
    index_map = MOI.copy_to(cache, src)
    A = convert(
        SparseArrays.SparseMatrixCSC{Float64,Int},
        cache.constraints.coefficients,
    )
    b = cache.constraints.constants
    c = zeros(size(A, 2))
    offset = cache.objective.scalar_affine.constant
    for term in cache.objective.scalar_affine.terms
        c[term.variable.value] += term.coefficient
    end
    if cache.objective.sense == MOI.MAX_SENSE
        c *= -1
    end
    xl, xu = cache.variables.lower, cache.variables.upper
    dest.termination_status, x_primal = solve_pdhg(A, b, c, xl, xu)
    for x in MOI.get(src, MOI.ListOfVariableIndices())
        dest.x_primal[x] = x_primal[index_map[x].value]
    end
    return index_map, false
end

```

If you’re wanting to work on this, I think you’ll need to dig and learn more of the MathOptInterface internals, including the definitions of what the various sets and functions are. Unfortunately, much of the `@product_of_sets` stuff is not well documented, because it isn’t meant for widespread public consumption. It’s really just an internal helper for solver wrappers, with the assumption that people using it are familiar with MOI. The best place to look are solver wrappers like Cbc.jl, Clp.jl, and SCS.jl.

For example,

```plaintext
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex{Float64}},
    ::Type{MOI.Interval},
)
    return true
end

```

needs to be

```plaintext
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{MOI.Interval{Float64}},
)
    return true
end

```

`VariableIndex` does not have a type parameter, `Interval` does.

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 15, 2023, 11:04pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/24 "2023-08-15T23:04:06Z")

</div>

@odow

Here is the template that I ran

I printed the time stamp right before `optimize!` and the timestamp at the first line of `optimize!` and it took 31 Seconds. For a larger problem, it took 7 minutes. At the same time, Mosek took 7 minutes to solve the large problem.

```julia
22:51:16 Calling optimize!...
22:51:47 GENERIC Solver Caching ...

```

[Data file](https://drive.google.com/file/d/105jCmzAnrcsbWVPtHPxKYV1akITpw7Q0/view?usp=sharing)

```julia
module GENERIC

using SparseArrays
using Dates
using LinearAlgebra
import MathOptInterface as MOI

MOI.Utilities.@product_of_sets(RHS, MOI.Nonpositives)

const OptimizerCache = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            Int,
            MOI.Utilities.OneBasedIndexing,
        },
        Vector{Float64},
        RHS{Float64},
    },
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    x_primal::Dict{MOI.VariableIndex, Float64}
    termination_status::MOI.TerminationStatusCode

    function Optimizer(; kwargs...)
        return new(Dict{MOI.VariableIndex, Float64}(), MOI.OPTIMIZE_NOT_CALLED)
    end
end

MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.is_empty(model::Optimizer)
    return isempty(model.x_primal) &&
        model.termination_status == MOI.OPTIMIZE_NOT_CALLED
end

function MOI.empty!(model::Optimizer)
    empty!(model.x_primal)
    model.termination_status = MOI.OPTIMIZE_NOT_CALLED
    return
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{MOI.Nonpositives},
)
    return true
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}},
)
    return true
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver Caching ...")
    cache = OptimizerCache()
    index_map = MOI.copy_to(cache, src)
    @assert all(==(-Inf), cache.variables.lower)
    @assert all(==(Inf), cache.variables.upper)

    println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver A and b ...")
    A = convert(
        SparseArrays.SparseMatrixCSC{Float64,Int},
        cache.constraints.coefficients,
    )
    rows, cols = size(A)
    b = cache.constraints.constants * -1 # Flip sign because of Ax + b <= 0 formulation of MOI   
    c = zeros(cols)

    obj_type = MOI.get(src, MOI.ObjectiveFunctionType())
    if obj_type == MOI.ScalarAffineFunction{Float64}
        offset = cache.objective.scalar_affine.constant

        println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver c ...")
        for term in cache.objective.scalar_affine.terms
            c[term.variable.value] += term.coefficient
        end
        if cache.objective.sense == MOI.MAX_SENSE
            c *= -1
        end
        println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver ...")
    elseif obj_type == MOI.ScalarQuadraticFunction{Float64}
        offset = cache.objective.scalar_quadratic.constant

        for term in cache.objective.scalar_quadratic.affine_terms
            c[term.variable.value] += term.coefficient
        end
        if cache.objective.sense == MOI.MAX_SENSE
            c *= -1
        end

        Q = sparse(Symmetric(zeros(cols, cols)))
        for term in cache.objective.scalar_quadratic.quadratic_terms
            Q[term.variable_1.value, term.variable_2.value] += term.coefficient
        end    
        println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver ...")
    else
        throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction{obj_type}()))
    end
    return index_map, false
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, x::MOI.VariableIndex)
    return model.x_primal[x]
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return model.termination_status == MOI.OPTIMAL ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return "$(model.termination_status)"
end

MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.termination_status

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if model.termination_status == MOI.OPTIMAL
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

MOI.get(model::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION

function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
    setfield!(model.solver_settings, Symbol(param.name), value)
    return nothing
end

function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
    if in(Symbol(param.name), fieldnames(typeof(model.solver_settings)))
        return getfield(model.solver_settings, Symbol(param.name))
    end
    error("RawOptimizerAttribute with name $(param.name) is not set.")
end

MOI.get(::Optimizer, ::MOI.SolverName) = "GENERIC"

end

using JuMP, Dates
model = read_from_file("model.mps.bz2", format = MOI.FileFormats.FORMAT_REW)
set_optimizer(model, GENERIC.Optimizer)
println(Dates.format(now(), "HH:MM:SS"), " Calling optimize!...")
optimize!(model)

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 15, 2023, 11:28pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/25 "2023-08-15T23:28:10Z")

</div>

This is a very dense problem with `15_408_073` non-zeros, and every constraint is `ScalarAffineFunction-in-LessThan`, so you’re forcing JuMP to bridge to `VectorAffineFunction-in-Nonpositives`.

Here’s a better way to do it, which takes only one second (although reading from file take a while):

```plaintext
julia> module GENERIC

       using SparseArrays
       using Dates
       using LinearAlgebra
       import MathOptInterface as MOI

       MOI.Utilities.@product_of_sets(RHS, MOI.LessThan{T})

       const OptimizerCache = MOI.Utilities.GenericModel{
           Float64,
           MOI.Utilities.ObjectiveContainer{Float64},
           MOI.Utilities.VariablesContainer{Float64},
           MOI.Utilities.MatrixOfConstraints{
               Float64,
               MOI.Utilities.MutableSparseMatrixCSC{
                   Float64,
                   Int,
                   MOI.Utilities.OneBasedIndexing,
               },
               MOI.Utilities.Hyperrectangle{Float64},
               RHS{Float64},
           },
       }

       mutable struct Optimizer <: MOI.AbstractOptimizer
           x_primal::Dict{MOI.VariableIndex, Float64}
           termination_status::MOI.TerminationStatusCode

           function Optimizer(; kwargs...)
               return new(Dict{MOI.VariableIndex, Float64}(), MOI.OPTIMIZE_NOT_CALLED)
           end
       end

       MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

       function MOI.is_empty(model::Optimizer)
           return isempty(model.x_primal) &&
               model.termination_status == MOI.OPTIMIZE_NOT_CALLED
       end

       function MOI.empty!(model::Optimizer)
           empty!(model.x_primal)
           model.termination_status = MOI.OPTIMIZE_NOT_CALLED
           return
       end

       function MOI.supports_constraint(
           ::Optimizer,
           ::Type{MOI.ScalarAffineFunction{Float64}},
           ::Type{MOI.LessThan{Float64}},
       )
           return true
       end

       function MOI.supports(
           ::Optimizer,
           ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
       )
           return true
       end

       function MOI.supports(
           ::Optimizer,
           ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}},
       )
           return true
       end

       function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
           println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver Caching ...")
           cache = OptimizerCache()
           index_map = MOI.copy_to(cache, src)
           println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver A and b ...")
           A = convert(
               SparseArrays.SparseMatrixCSC{Float64,Int},
               cache.constraints.coefficients,
           )
           rows, cols = size(A)
           b = cache.constraints.constants.upper
           c = zeros(cols)
           obj_type = MOI.get(src, MOI.ObjectiveFunctionType())
           @assert obj_type == MOI.ScalarAffineFunction{Float64}
           offset = cache.objective.scalar_affine.constant

           println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver c ...")
           for term in cache.objective.scalar_affine.terms
               c[term.variable.value] += term.coefficient
           end
           if cache.objective.sense == MOI.MAX_SENSE
               c *= -1
           end
           println(Dates.format(now(), "HH:MM:SS"), " GENERIC Solver ...")
           return index_map, false
       end

       function MOI.get(model::Optimizer, ::MOI.VariablePrimal, x::MOI.VariableIndex)
           return model.x_primal[x]
       end

       function MOI.get(model::Optimizer, ::MOI.ResultCount)
           return model.termination_status == MOI.OPTIMAL ? 1 : 0
       end

       function MOI.get(model::Optimizer, ::MOI.RawStatusString)
           return "$(model.termination_status)"
       end

       MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.termination_status

       function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
           if model.termination_status == MOI.OPTIMAL
               return MOI.FEASIBLE_POINT
           else
               return MOI.NO_SOLUTION
           end
       end

       MOI.get(model::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION

       MOI.get(::Optimizer, ::MOI.SolverName) = "GENERIC"

       end
WARNING: replacing module GENERIC.
Main.GENERIC

julia> using JuMP, Dates

julia> @time model = read_from_file("model.mps.bz2", format = MOI.FileFormats.FORMAT_REW)
 50.414692 seconds (140.04 M allocations: 11.260 GiB, 3.97% gc time)
A JuMP Model
Minimization problem with:
Variables: 17682
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 17706 constraints
`VariableRef`-in-`MathOptInterface.EqualTo{Float64}`: 1 constraint
`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 10000 constraints
`VariableRef`-in-`MathOptInterface.Interval{Float64}`: 7680 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> set_optimizer(model, GENERIC.Optimizer; add_bridges = false)

julia> println(Dates.format(now(), "HH:MM:SS"), " Calling optimize!...")
11:26:21 Calling optimize!...

julia> optimize!(model)
11:26:21 GENERIC Solver Caching ...
11:26:22 GENERIC Solver A and b ...
11:26:22 GENERIC Solver c ...
11:26:22 GENERIC Solver ...

```

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 15, 2023, 11:53pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/26 "2023-08-15T23:53:42Z")

</div>

Thanks, it definitely helped.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 16, 2023, 12:16am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/27 "2023-08-16T00:16:29Z")

</div>

What’s your goal for this? Are you just trying to write a parser to convert MPS files into the data matrices?

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 16, 2023, 12:29am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/28 "2023-08-16T00:29:32Z")

</div>

No, integrating an internal new LP/QP solver on GPU. I transferred the problem into an MPS so I could share it.

When I connect directly now, I get

```julia
ERROR: LoadError: Constraints of type MathOptInterface.VariableIndex-in-MathOptInterface.GreaterThan{Float64} are not supported by the solver.

```

I guess, as you said, I need to read MOI thoroughly.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 16, 2023, 12:34am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/29 "2023-08-16T00:34:09Z")

</div>

I cannot reproduce your error:

```julia
julia> using JuMP

julia> model = Model(GENERIC.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: GENERIC

julia> @variable(model, x >= 0)
x

julia> @objective(model, Min, 2x)
2 x

julia> optimize!(model)
12:32:53 GENERIC Solver Caching ...
12:32:53 GENERIC Solver A and b ...
12:32:53 GENERIC Solver c ...
12:32:53 GENERIC Solver ...

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 16, 2023, 12:37am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/30 "2023-08-16T00:37:36Z")

</div>

Connecting a new solver to JuMP is quite a lot of work:

[https://jump.dev/MathOptInterface.jl/stable/tutorials/implementing/](https://jump.dev/MathOptInterface.jl/stable/tutorials/implementing/)

The easiest approach is to find a similar solver already implemented. The code above was mainly intended as a simple interface for a toy solver. You probably need quite a bit more for a proper interface. (The “simple” Clp wrapper is 800 lines: [https://github.com/jump-dev/Clp.jl/blob/master/src/MOI\_wrapper/MOI\_wrapper.jl](https://github.com/jump-dev/Clp.jl/blob/master/src/MOI_wrapper/MOI_wrapper.jl))

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 16, 2023, 12:38am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/31 "2023-08-16T00:38:32Z")

</div>

When I exported the problem to MPS, it may represent l≤x≤b as MOI.LessThan (x≤b, -x≤-l). My code is integrated with Python, so mimicking the code here will be very difficult.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 16, 2023, 12:43am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/32 "2023-08-16T00:43:59Z")

</div>

> When I exported the problem to MPS

Not sure I follow. JuMP doesn’t reformulate `l <= x <= b` to `-x <= -l` when writing to MPS. And I don’t know if I understand the Python integration part.

Is the solver written in Julia? Why are you trying to integrate with JuMP if you’re calling it from Python? Is the intention for it to consume MPS files? Or do you want to solve models written directly in JuMP?

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 16, 2023, 12:48am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/33 "2023-08-16T00:48:36Z")

</div>

Solver is in C++, and I used CxxWrap.jl to create .so to link it with Julia. Python code provides the forecasting and calls the Julia JuMP via PyJulia.

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 16, 2023, 12:49am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/34 "2023-08-16T00:49:55Z")

</div>

Anyhow I need to figure out how to add this for now `MathOptInterface.VariableIndex-in-MathOptInterface.GreaterThan{Float64}`

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 16, 2023, 2:40am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/35 "2023-08-16T02:40:35Z")

</div>

Sometimes JuMP is too clever for its own good. I forgot the `supports_constraint` for variables, so it was bridging to `-x <= -l`.

```julia
module GENERIC

import Dates
import MathOptInterface as MOI
import SparseArrays

MOI.Utilities.@product_of_sets(ConstraintSets, MOI.LessThan{T})

const OptimizerCache = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            Int,
            MOI.Utilities.OneBasedIndexing,
        },
        MOI.Utilities.Hyperrectangle{Float64},
        ConstraintSets{Float64},
    },
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    x_primal::Dict{MOI.VariableIndex, Float64}
    termination_status::MOI.TerminationStatusCode

    function Optimizer(; kwargs...)
        return new(Dict{MOI.VariableIndex, Float64}(), MOI.OPTIMIZE_NOT_CALLED)
    end
end

MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.is_empty(model::Optimizer)
    return isempty(model.x_primal) &&
        model.termination_status == MOI.OPTIMIZE_NOT_CALLED
end

function MOI.empty!(model::Optimizer)
    empty!(model.x_primal)
    model.termination_status = MOI.OPTIMIZE_NOT_CALLED
    return
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.EqualTo{Float64},
            MOI.Interval{Float64},
        }
    },
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.LessThan{Float64}},
)
    return true
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{
        <:Union{
            MOI.ScalarAffineFunction{Float64},
            MOI.ScalarQuadraticFunction{Float64},
        },
    },
)
    return true
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    println(Dates.format(Dates.now(), "HH:MM:SS"), " GENERIC Solver Caching ...")
    cache = OptimizerCache()
    index_map = MOI.copy_to(cache, src)
    println(Dates.format(Dates.now(), "HH:MM:SS"), " GENERIC Solver A and b ...")
    A = convert(
        SparseArrays.SparseMatrixCSC{Float64,Int},
        cache.constraints.coefficients,
    )
    rows, cols = size(A)
    c = zeros(cols)
    obj_type = MOI.get(src, MOI.ObjectiveFunctionType())
    @assert obj_type == MOI.ScalarAffineFunction{Float64}
    println(Dates.format(Dates.now(), "HH:MM:SS"), " GENERIC Solver c ...")
    for term in cache.objective.scalar_affine.terms
        c[term.variable.value] += term.coefficient
    end
    if cache.objective.sense == MOI.MAX_SENSE
        c *= -1
    end
    println(Dates.format(Dates.now(), "HH:MM:SS"), " GENERIC Solver ...")
    # solve(
    # A, 
    # cache.constraints.constants.upper,
    # c, 
    # cache.objective.scalar_affine.constant,
    # cache.variables.lower,
    # cache.variables.upper,
    # )
    return index_map, false
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, x::MOI.VariableIndex)
    return model.x_primal[x]
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return model.termination_status == MOI.OPTIMAL ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return "$(model.termination_status)"
end

MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.termination_status

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if model.termination_status == MOI.OPTIMAL
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

MOI.get(model::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION

MOI.get(::Optimizer, ::MOI.SolverName) = "GENERIC"

end

using JuMP
model = Model(GENERIC.Optimizer; add_bridges = false)
@variable(model, x >= 0)
@objective(model, Min, 2x)
optimize!(model)

```

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 16, 2023, 3:27am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/36 "2023-08-16T03:27:07Z")

</div>

Thank you. Now I can solve it.

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 16, 2023, 6:21am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/37 "2023-08-16T06:21:19Z")

</div>

I also Added, I hope it is correct.

```julia
    finite_upper = .!isinf.(cache.variables.upper)
    finite_lower = .!isinf.(cache.variables.lower)
    
    len_upper = sum(finite_upper)
    len_lower = sum(finite_lower)
    num_columns = cache.constraints.coefficients.n
    
    A = vcat(
        convert(SparseArrays.SparseMatrixCSC{Float64, Int}, cache.constraints.coefficients),
        sparse(1:len_upper, findall(finite_upper), 1.0, len_upper, num_columns),
        sparse(1:len_lower, findall(finite_lower), -1.0, len_lower, num_columns)
    )

    b = vcat(
        cache.constraints.constants.upper,
        cache.variables.upper[finite_upper],
        cache.variables.lower[finite_lower]
    )

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 16, 2023, 7:39am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/38 "2023-08-16T07:39:40Z")

</div>

So your solver _doesn’t_ support variable bounds? Why not just let JuMP reformulate things?

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 16, 2023, 3:51pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/39 "2023-08-16T15:51:28Z")

</div>

Using bridge? Because it takes too long.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 16, 2023, 8:14pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/40 "2023-08-16T20:14:48Z")

</div>

You could let JuMP reformulate, and pay the extra few seconds, or you could reformulate your model to be closer to the solver.

Instead of

```plaintext
model = Model(GENERIC.Optimizer)
@variable(model, x >= 0)

```

do

```plaintext
model = Model(GENERIC.Optimizer)
@variable(model, x)
@constraint(model, x >= 0)

```

See [Ellipsoid approximation · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/conic/ellipse_approx/#Alternative-formulations) for some explanation.

[Previous page](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694.md?page=1)