# Load JuMP
using JuMP
using MathOptInterface
# Load solver package
using GLPK
# shortcuts
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

# Use a hyperplane representation of the polyhedron
# i.e. P = { x | Ax ≤ b }
n = 4
# Store LHS as vector-of-vectors
A = Vector{Float64}[
    [ 2, 1], [ 2,-1],
    [-1, 2], [-1,-2]]
b = ones(n)

# Build JuMP model
m = Model(optimizer = GLPK.GLPKOptimizerLP())
@variable(m, r)
@variable(m, x[1:2])
@objective(m, Max, 1r)
@constraint(m, P[i=1:4], sum(A[i][j]*x[j] for j in eachindex(x)) + norm(A[i])*r <= b[i])
# solve problem
JuMP.optimize(m)

# Where is the center?
print(JuMP.resultvalue.(x))

JuMP.resultvalue(r)

# Use Gadfly to display the solution
using Gadfly
Gadfly.set_default_plot_size(8cm, 8cm)
# Plot lines over [-1.5, 1.5]
domain = linspace(-1, +1)
# Plot circle across all angles
θ = linspace(0,2π)
plot(
# a_1 x_1 + a_2 x_2 = b
# --> x_2 = (b - a_1 x_1)/a_2
[layer(x=domain,
       y=(b[i]-A[i][1]*domain)/A[i][2],
       Geom.line,
       Theme(line_width=2px,
             default_color=colorant"blue")) for i in 1:4]...,
# x_1' = x_1 + rθ
# x_2' = x_2 + rθ
layer(x=JuMP.resultvalue(x[1]) + JuMP.resultvalue(r)*cos.(θ),
      y=JuMP.resultvalue(x[2]) + JuMP.resultvalue(r)*sin.(θ),
      Geom.path,
      Theme(line_width=5px,
            default_color=colorant"red")),
Coord.Cartesian(ymin=-1,ymax=+1)
)