Skip to content

Instantly share code, notes, and snippets.

@jonocarroll

jonocarroll/mip_julia.jl

Last active October 22, 2019 12:28
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jonocarroll/f00036b9c00bebf3923d64f1a8ffca91 to your computer and use it in GitHub Desktop.
Save jonocarroll/f00036b9c00bebf3923d64f1a8ffca91 to your computer and use it in GitHub Desktop.
Download ZIP
Graph coloring comparison between Julia and R
Raw
mip_julia.jl
using JuMP, GLPK
function graph_coloring(;fast=true, benchmark=false)
m = Model(with_optimizer(GLPK.Optimizer, msg_lev = 3))
@variable(m, max_color, Int)
@variable(m, c[1:49,1:4], Bin)
@objective(m, Min, max_color)
for s=1:49
@constraint(m, sum(c[s,1:4]) == 1)
end
if fast
for s=1:49
@constraint(m, sum(c[s,i]*i for i=1:4) <= max_color)
end
else
for s=1:49, i=1:4
@constraint(m, max_color >= i*c[s,i])
end
end
for i=1:4
@constraint(m, c[1,i] + c[11,i] <= 1)
@constraint(m, c[1,i] + c[8,i] <= 1)
@constraint(m, c[10,i] + c[13,i] <= 1)
@constraint(m, c[10,i] + c[7,i] <= 1)
@constraint(m, c[11,i] + c[22,i] <= 1)
@constraint(m, c[11,i] + c[24,i] <= 1)
@constraint(m, c[11,i] + c[8,i] <= 1)
@constraint(m, c[12,i] + c[14,i] <= 1)
@constraint(m, c[12,i] + c[9,i] <= 1)
@constraint(m, c[13,i] + c[26,i] <= 1)
@constraint(m, c[13,i] + c[35,i] <= 1)
@constraint(m, c[13,i] + c[7,i] <= 1)
@constraint(m, c[14,i] + c[16,i] <= 1)
@constraint(m, c[14,i] + c[18,i] <= 1)
@constraint(m, c[14,i] + c[19,i] <= 1)
@constraint(m, c[15,i] + c[13,i] <= 1)
@constraint(m, c[15,i] + c[33,i] <= 1)
@constraint(m, c[15,i] + c[35,i] <= 1)
@constraint(m, c[16,i] + c[17,i] <= 1)
@constraint(m, c[16,i] + c[18,i] <= 1)
@constraint(m, c[16,i] + c[20,i] <= 1)
@constraint(m, c[17,i] + c[20,i] <= 1)
@constraint(m, c[17,i] + c[28,i] <= 1)
@constraint(m, c[17,i] + c[29,i] <= 1)
@constraint(m, c[17,i] + c[30,i] <= 1)
@constraint(m, c[18,i] + c[19,i] <= 1)
@constraint(m, c[2,i] + c[4,i] <= 1)
@constraint(m, c[2,i] + c[5,i] <= 1)
@constraint(m, c[2,i] + c[6,i] <= 1)
@constraint(m, c[20,i] + c[28,i] <= 1)
@constraint(m, c[21,i] + c[25,i] <= 1)
@constraint(m, c[22,i] + c[23,i] <= 1)
@constraint(m, c[22,i] + c[36,i] <= 1)
@constraint(m, c[22,i] + c[8,i] <= 1)
@constraint(m, c[23,i] + c[31,i] <= 1)
@constraint(m, c[23,i] + c[36,i] <= 1)
@constraint(m, c[23,i] + c[41,i] <= 1)
@constraint(m, c[23,i] + c[6,i] <= 1)
@constraint(m, c[24,i] + c[22,i] <= 1)
@constraint(m, c[24,i] + c[36,i] <= 1)
@constraint(m, c[25,i] + c[17,i] <= 1)
@constraint(m, c[25,i] + c[29,i] <= 1)
@constraint(m, c[26,i] + c[21,i] <= 1)
@constraint(m, c[26,i] + c[32,i] <= 1)
@constraint(m, c[27,i] + c[34,i] <= 1)
@constraint(m, c[29,i] + c[34,i] <= 1)
@constraint(m, c[3,i] + c[12,i] <= 1)
@constraint(m, c[30,i] + c[27,i] <= 1)
@constraint(m, c[30,i] + c[29,i] <= 1)
@constraint(m, c[30,i] + c[34,i] <= 1)
@constraint(m, c[31,i] + c[15,i] <= 1)
@constraint(m, c[31,i] + c[33,i] <= 1)
@constraint(m, c[31,i] + c[37,i] <= 1)
@constraint(m, c[31,i] + c[41,i] <= 1)
@constraint(m, c[32,i] + c[21,i] <= 1)
@constraint(m, c[32,i] + c[25,i] <= 1)
@constraint(m, c[32,i] + c[29,i] <= 1)
@constraint(m, c[32,i] + c[34,i] <= 1)
@constraint(m, c[32,i] + c[39,i] <= 1)
@constraint(m, c[33,i] + c[35,i] <= 1)
@constraint(m, c[33,i] + c[37,i] <= 1)
@constraint(m, c[34,i] + c[38,i] <= 1)
@constraint(m, c[35,i] + c[26,i] <= 1)
@constraint(m, c[35,i] + c[32,i] <= 1)
@constraint(m, c[35,i] + c[39,i] <= 1)
@constraint(m, c[35,i] + c[46,i] <= 1)
@constraint(m, c[36,i] + c[31,i] <= 1)
@constraint(m, c[36,i] + c[41,i] <= 1)
@constraint(m, c[37,i] + c[40,i] <= 1)
@constraint(m, c[37,i] + c[46,i] <= 1)
@constraint(m, c[38,i] + c[45,i] <= 1)
@constraint(m, c[39,i] + c[34,i] <= 1)
@constraint(m, c[39,i] + c[38,i] <= 1)
@constraint(m, c[39,i] + c[42,i] <= 1)
@constraint(m, c[39,i] + c[43,i] <= 1)
@constraint(m, c[39,i] + c[44,i] <= 1)
@constraint(m, c[4,i] + c[10,i] <= 1)
@constraint(m, c[4,i] + c[5,i] <= 1)
@constraint(m, c[40,i] + c[46,i] <= 1)
@constraint(m, c[40,i] + c[47,i] <= 1)
@constraint(m, c[41,i] + c[37,i] <= 1)
@constraint(m, c[41,i] + c[40,i] <= 1)
@constraint(m, c[42,i] + c[44,i] <= 1)
@constraint(m, c[42,i] + c[48,i] <= 1)
@constraint(m, c[43,i] + c[42,i] <= 1)
@constraint(m, c[44,i] + c[48,i] <= 1)
@constraint(m, c[45,i] + c[44,i] <= 1)
@constraint(m, c[46,i] + c[39,i] <= 1)
@constraint(m, c[46,i] + c[43,i] <= 1)
@constraint(m, c[46,i] + c[47,i] <= 1)
@constraint(m, c[47,i] + c[43,i] <= 1)
@constraint(m, c[49,i] + c[21,i] <= 1)
@constraint(m, c[49,i] + c[25,i] <= 1)
@constraint(m, c[5,i] + c[10,i] <= 1)
@constraint(m, c[5,i] + c[13,i] <= 1)
@constraint(m, c[5,i] + c[15,i] <= 1)
@constraint(m, c[6,i] + c[15,i] <= 1)
@constraint(m, c[6,i] + c[31,i] <= 1)
@constraint(m, c[6,i] + c[5,i] <= 1)
@constraint(m, c[7,i] + c[26,i] <= 1)
@constraint(m, c[7,i] + c[49,i] <= 1)
@constraint(m, c[8,i] + c[2,i] <= 1)
@constraint(m, c[8,i] + c[23,i] <= 1)
@constraint(m, c[8,i] + c[6,i] <= 1)
@constraint(m, c[9,i] + c[14,i] <= 1)
@constraint(m, c[9,i] + c[16,i] <= 1)
end
optimize!(m)
if !benchmark
println(JuMP.objective_value(m))
println(JuMP.termination_status(m))
end
end
using BenchmarkTools
@benchmark graph_coloring(fast=true,benchmark=false)
# GLPK Simplex Optimizer, v4.64
# 522 rows, 197 columns, 1289 non-zeros
# 0: obj = 0.000000000e+00 inf = 4.900e+01 (49)
# 117: obj = 1.500000000e+00 inf = 0.000e+00 (0)
# * 121: obj = 1.500000000e+00 inf = 0.000e+00 (0)
# OPTIMAL LP SOLUTION FOUND
# GLPK Integer Optimizer, v4.64
# 522 rows, 197 columns, 1289 non-zeros
# 197 integer variables, 196 of which are binary
# Integer optimization begins...
# + 121: mip = not found yet >= -inf (1; 0)
# + 408: >>>>> 4.000000000e+00 >= 2.000000000e+00 50.0% (31; 0)
# + 7008: mip = 4.000000000e+00 >= tree is empty 0.0% (0; 709)
# INTEGER OPTIMAL SOLUTION FOUND
# 4.0
# OPTIMAL
# BenchmarkTools.Trial:
# memory estimate: 2.43 MiB
# allocs estimate: 43346
# --------------
# minimum time: 621.202 ms (0.00% GC)
# median time: 644.946 ms (0.00% GC)
# mean time: 645.873 ms (0.00% GC)
# maximum time: 669.376 ms (0.00% GC)
# --------------
# samples: 8
# evals/sample: 1
Raw
mip_r.R
library(rmpk)
library(ROI.plugin.glpk)
# number of states to color
n <- 49
model <- MIPModel(ROI_solver("glpk", control = list(presolve = FALSE, verbose = TRUE)))
# @variable(m, max_color, Int)
model$add_variable(max_colors, type = "integer", lb = 1, ub = 10)
# @objective(m, Min, max_color)
model$set_objective(max_colors, sense = "min")
# @variable(m, c[1:49,1:4], Bin)
model$add_variable(x[i, k], type = "binary", i = 1:n, k = 1:4)
# for s=1:49
# @constraint(m, sum(c[s,1:4]) == 1)
# end
for (s in 1:n) {
model$add_constraint(sum_expr(x[s, k], k = 1:4) == 1)
}
# for s=1:49
# @constraint(m, sum(c[s,i]*i for i=1:4) <= max_color)
# end
for (s in 1:n) {
model$add_constraint(sum_expr(k * x[s, k], k = 1:4) <= max_colors)
}
# no adjacent nodes have the same color
for (k in 1:4) {
model$add_constraint(x[1,k] + x[11,k] <= 1)
model$add_constraint(x[1,k] + x[8,k] <= 1)
model$add_constraint(x[11,k] + x[8,k] <= 1)
model$add_constraint(x[11,k] + x[22,k] <= 1)
model$add_constraint(x[11,k] + x[24,k] <= 1)
model$add_constraint(x[24,k] + x[22,k] <= 1)
model$add_constraint(x[24,k] + x[36,k] <= 1)
model$add_constraint(x[22,k] + x[8,k] <= 1)
model$add_constraint(x[22,k] + x[23,k] <= 1)
model$add_constraint(x[22,k] + x[36,k] <= 1)
model$add_constraint(x[8,k] + x[2,k] <= 1)
model$add_constraint(x[8,k] + x[6,k] <= 1)
model$add_constraint(x[8,k] + x[23,k] <= 1)
model$add_constraint(x[23,k] + x[6,k] <= 1)
model$add_constraint(x[23,k] + x[31,k] <= 1)
model$add_constraint(x[23,k] + x[41,k] <= 1)
model$add_constraint(x[23,k] + x[36,k] <= 1)
model$add_constraint(x[36,k] + x[31,k] <= 1)
model$add_constraint(x[36,k] + x[41,k] <= 1)
model$add_constraint(x[2,k] + x[4,k] <= 1)
model$add_constraint(x[2,k] + x[5,k] <= 1)
model$add_constraint(x[2,k] + x[6,k] <= 1)
model$add_constraint(x[6,k] + x[5,k] <= 1)
model$add_constraint(x[6,k] + x[15,k] <= 1)
model$add_constraint(x[6,k] + x[31,k] <= 1)
model$add_constraint(x[31,k] + x[15,k] <= 1)
model$add_constraint(x[31,k] + x[33,k] <= 1)
model$add_constraint(x[31,k] + x[37,k] <= 1)
model$add_constraint(x[31,k] + x[41,k] <= 1)
model$add_constraint(x[41,k] + x[37,k] <= 1)
model$add_constraint(x[41,k] + x[40,k] <= 1)
model$add_constraint(x[4,k] + x[10,k] <= 1)
model$add_constraint(x[4,k] + x[5,k] <= 1)
model$add_constraint(x[5,k] + x[10,k] <= 1)
model$add_constraint(x[5,k] + x[13,k] <= 1)
model$add_constraint(x[5,k] + x[15,k] <= 1)
model$add_constraint(x[15,k] + x[13,k] <= 1)
model$add_constraint(x[15,k] + x[35,k] <= 1)
model$add_constraint(x[15,k] + x[33,k] <= 1)
model$add_constraint(x[33,k] + x[35,k] <= 1)
model$add_constraint(x[33,k] + x[37,k] <= 1)
model$add_constraint(x[37,k] + x[46,k] <= 1)
model$add_constraint(x[37,k] + x[40,k] <= 1)
model$add_constraint(x[40,k] + x[46,k] <= 1)
model$add_constraint(x[40,k] + x[47,k] <= 1)
model$add_constraint(x[10,k] + x[7,k] <= 1)
model$add_constraint(x[10,k] + x[13,k] <= 1)
model$add_constraint(x[13,k] + x[7,k] <= 1)
model$add_constraint(x[13,k] + x[26,k] <= 1)
model$add_constraint(x[13,k] + x[35,k] <= 1)
model$add_constraint(x[35,k] + x[26,k] <= 1)
model$add_constraint(x[35,k] + x[32,k] <= 1)
model$add_constraint(x[35,k] + x[39,k] <= 1)
model$add_constraint(x[35,k] + x[46,k] <= 1)
model$add_constraint(x[46,k] + x[39,k] <= 1)
model$add_constraint(x[46,k] + x[43,k] <= 1)
model$add_constraint(x[46,k] + x[47,k] <= 1)
model$add_constraint(x[47,k] + x[43,k] <= 1)
model$add_constraint(x[7,k] + x[26,k] <= 1)
model$add_constraint(x[7,k] + x[49,k] <= 1)
model$add_constraint(x[26,k] + x[21,k] <= 1)
model$add_constraint(x[26,k] + x[32,k] <= 1)
model$add_constraint(x[32,k] + x[21,k] <= 1)
model$add_constraint(x[32,k] + x[25,k] <= 1)
model$add_constraint(x[32,k] + x[29,k] <= 1)
model$add_constraint(x[32,k] + x[34,k] <= 1)
model$add_constraint(x[32,k] + x[39,k] <= 1)
model$add_constraint(x[39,k] + x[34,k] <= 1)
model$add_constraint(x[39,k] + x[38,k] <= 1)
model$add_constraint(x[39,k] + x[44,k] <= 1)
model$add_constraint(x[39,k] + x[42,k] <= 1)
model$add_constraint(x[39,k] + x[43,k] <= 1)
model$add_constraint(x[43,k] + x[42,k] <= 1)
model$add_constraint(x[49,k] + x[21,k] <= 1)
model$add_constraint(x[49,k] + x[25,k] <= 1)
model$add_constraint(x[21,k] + x[25,k] <= 1)
model$add_constraint(x[42,k] + x[44,k] <= 1)
model$add_constraint(x[42,k] + x[48,k] <= 1)
model$add_constraint(x[25,k] + x[17,k] <= 1)
model$add_constraint(x[25,k] + x[29,k] <= 1)
model$add_constraint(x[3,k] + x[12,k] <= 1)
model$add_constraint(x[12,k] + x[9,k] <= 1)
model$add_constraint(x[12,k] + x[14,k] <= 1)
model$add_constraint(x[9,k] + x[14,k] <= 1)
model$add_constraint(x[9,k] + x[16,k] <= 1)
model$add_constraint(x[14,k] + x[16,k] <= 1)
model$add_constraint(x[14,k] + x[19,k] <= 1)
model$add_constraint(x[14,k] + x[18,k] <= 1)
model$add_constraint(x[18,k] + x[19,k] <= 1)
model$add_constraint(x[16,k] + x[18,k] <= 1)
model$add_constraint(x[16,k] + x[17,k] <= 1)
model$add_constraint(x[16,k] + x[20,k] <= 1)
model$add_constraint(x[17,k] + x[20,k] <= 1)
model$add_constraint(x[17,k] + x[28,k] <= 1)
model$add_constraint(x[17,k] + x[30,k] <= 1)
model$add_constraint(x[17,k] + x[29,k] <= 1)
model$add_constraint(x[20,k] + x[28,k] <= 1)
model$add_constraint(x[30,k] + x[27,k] <= 1)
model$add_constraint(x[30,k] + x[29,k] <= 1)
model$add_constraint(x[30,k] + x[34,k] <= 1)
model$add_constraint(x[27,k] + x[34,k] <= 1)
model$add_constraint(x[29,k] + x[34,k] <= 1)
model$add_constraint(x[34,k] + x[38,k] <= 1)
model$add_constraint(x[38,k] + x[45,k] <= 1)
model$add_constraint(x[45,k] + x[44,k] <= 1)
model$add_constraint(x[44,k] + x[48,k] <= 1)
}
model
# MIP Model:
# Variables: 197
# Constraints: 522
microbenchmark::microbenchmark(model$optimize())
## <SOLVER MSG> ----
## GLPK Simplex Optimizer, v4.65
## 522 rows, 197 columns, 1289 non-zeros
## 0: obj = 1.000000000e+00 inf = 4.900e+01 (49)
## 115: obj = 1.500000000e+00 inf = 0.000e+00 (0)
## * 116: obj = 1.500000000e+00 inf = 0.000e+00 (0)
## OPTIMAL LP SOLUTION FOUND
## GLPK Integer Optimizer, v4.65
## 522 rows, 197 columns, 1289 non-zeros
## 197 integer variables, 196 of which are binary
## Integer optimization begins...
## Long-step dual simplex will be used
## + 116: mip = not found yet >= -inf (1; 0)
## + 388: >>>>> 4.000000000e+00 >= 2.000000000e+00 50.0% (24; 0)
## + 3033: mip = 4.000000000e+00 >= tree is empty 0.0% (0; 501)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----
## Unit: milliseconds
## expr min lq mean median uq max neval
## model$optimize() 325.6906 341.3275 350.1407 347.463 356.3246 402.5517 100
model$objective_value()
## 4
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment