Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Demo code for Stéphan Laurent blog post comment
import GLPK # NOTE: it's better to use 'import' with GLPK,
# rather than 'using'
μ = [1/7, 2/7, 4/7]
ν = [1/4, 1/4, 1/2]
n = length(μ)
D = 1 .- eye(n) # NOTE: starting from julia v0.3, you won't be
# allowed to use '-' here, you need to
# use '.-'
# alternative definition, only works in julia v0.3:
# D = ones(n,n) - I
M1 = zeros(n, n^2)
for i = 1:n
M1[i, (1:n)+n*(i-1)] = 1
end
M2 = repmat(eye(n), 1, n)
A = vcat(M1, M2)
lp = GLPK.Prob()
GLPK.set_prob_name(lp, "kanto")
GLPK.set_obj_dir(lp, GLPK.MIN)
GLPK.add_rows(lp, 2n)
for i = 1:n
GLPK.set_row_bnds(lp, i, GLPK.FX, μ[i], μ[i])
GLPK.set_row_bnds(lp, i+n, GLPK.FX, ν[i], ν[i])
end
GLPK.add_cols(lp, n^2)
k = 0
for i = 1:n, j = 1:n # NOTE: nested loops can be written like this
k += 1
GLPK.set_col_bnds(lp, k, GLPK.LO, 0.0, 0.0)
GLPK.set_obj_coef(lp, k, D[i,j])
end
GLPK.load_matrix(lp, sparse(A)) # NOTE: GLPK.load_matrix can work with sparse matrices
flag = GLPK.simplex(lp)
# flag = GLPK.exact(lp)
flag == 0 || error("GLPK failed")
result = GLPK.get_obj_val(lp)
using MathProgBase
# using GLPKMathProgInterface # NOTE: only needed to specify solver options
# e.g. method = :Exact
μ = [1/7, 2/7, 4/7]
ν = [1/4, 1/4, 1/2]
n = length(μ)
D = 1 .- eye(n) # NOTE: starting from julia v0.3, you won't be
# allowed to use '-' here, you need to
# use '.-'
# alternative definition, only works in julia v0.3:
# D = ones(n,n) - I
M1 = zeros(n, n^2)
for i = 1:n
M1[i, (1:n)+n*(i-1)] = 1
end
M2 = repmat(eye(n), 1, n)
A = vcat(M1, M2)
μν = vcat(μ, ν)
vD = vec(D)
sol = linprog(vD, A, '=', μν, 0.0, Inf)
#sol = linprog(vD, A, '=', μν, 0.0, Inf, GLPKSolverLP(method = :Exact))
sol.status == :Optimal || error("Optimization failed: $(sol.status)")
result = sol.objval
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment