diff --git a/LICENSE.md b/LICENSE.md
index a5fa124..3c472ad 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -1,7 +1,8 @@
 The MOIPaperBenchmark repository is licensed under the MIT "Expat" License:
 
-> Copyright (c) 2017: Benoît Legat, Oscar Dowson, Joaquim Dias Garcia,
-> Miles Lubin and contributors
+> Copyright (c) 2020: Benoît Legat, Oscar Dowson, Joaquim Dias Garcia,
+> and contributors
+> Copyright (c) 2020: Google LLC
 >
 > Permission is hereby granted, free of charge, to any person obtaining a copy
 > of this software and associated documentation files (the "Software"), to deal
diff --git a/Manifest.toml b/Manifest.toml
new file mode 100644
index 0000000..256903f
--- /dev/null
+++ b/Manifest.toml
@@ -0,0 +1,251 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[BenchmarkTools]]
+deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
+git-tree-sha1 = "9e62e66db34540a0c919d72172cc2f642ac71260"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "0.5.0"
+
+[[BinaryProvider]]
+deps = ["Libdl", "Logging", "SHA"]
+git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
+uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
+version = "0.5.10"
+
+[[Bzip2_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "3663bfffede2ef41358b6fc2e1d8a6d50b3c3904"
+uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
+version = "1.0.6+2"
+
+[[CEnum]]
+git-tree-sha1 = "1b77a77c3b28e0b3f413f7567c9bb8dd9bdccd14"
+uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
+version = "0.3.0"
+
+[[CodecBzip2]]
+deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
+git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
+uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
+version = "0.7.2"
+
+[[CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.0"
+
+[[CompilerSupportLibraries_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612"
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "0.3.3+0"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[GLPK]]
+deps = ["BinaryProvider", "CEnum", "GLPK_jll", "Libdl", "MathOptInterface"]
+git-tree-sha1 = "84afecdf40dd0102d798d24f63c5042773ad80f5"
+repo-rev = "master"
+repo-url = "https://github.com/jump-dev/GLPK.jl.git"
+uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+version = "0.14.0"
+
+[[GLPK_jll]]
+deps = ["GMP_jll", "Libdl", "Pkg"]
+git-tree-sha1 = "ccc855de74292e478d4278e3a6fdd8212f75e81e"
+uuid = "e8aa6df9-e6ca-548a-97ff-1f85fc5b8b98"
+version = "4.64.0+0"
+
+[[GMP_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "4dd9301d3a027c05ec403e756ee7a60e3c367e5d"
+uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
+version = "6.1.2+5"
+
+[[HTTP]]
+deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
+git-tree-sha1 = "eca61b35cdd8cd2fcc5eec1eda766424a995b02f"
+uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
+version = "0.8.16"
+
+[[IniFile]]
+deps = ["Test"]
+git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
+uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
+version = "0.5.0"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.0"
+
+[[JSONSchema]]
+deps = ["HTTP", "JSON", "ZipFile"]
+git-tree-sha1 = "832a4d327d9dafdae55a6ecae04f9997c83615cc"
+uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
+version = "0.3.0"
+
+[[LibGit2]]
+deps = ["Printf"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[MathOptInterface]]
+deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"]
+git-tree-sha1 = "cd2049c055c7d192a235670d50faa375361624ba"
+uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+version = "0.9.14"
+
+[[MathProgBase]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "9abbe463a1e9fc507f12a69e7f29346c2cdc472c"
+uuid = "fdba3010-5040-5b88-9595-932c9decdf73"
+version = "0.7.8"
+
+[[MbedTLS]]
+deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
+git-tree-sha1 = "426a6978b03a97ceb7ead77775a1da066343ec6e"
+uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
+version = "1.0.2"
+
+[[MbedTLS_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "a0cb0d489819fa7ea5f9fa84c7e7eba19d8073af"
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.16.6+1"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[MutableArithmetics]]
+deps = ["LinearAlgebra", "SparseArrays", "Test"]
+git-tree-sha1 = "6cf09794783b9de2e662c4e8b60d743021e338d0"
+uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+version = "0.2.10"
+
+[[OpenBLAS_jll]]
+deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
+git-tree-sha1 = "0c922fd9634e358622e333fc58de61f05a048492"
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.9+5"
+
+[[OrderedCollections]]
+git-tree-sha1 = "293b70ac1780f9584c89268a6e2a560d938a7065"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.3.0"
+
+[[Parsers]]
+deps = ["Dates", "Test"]
+git-tree-sha1 = "10134f2ee0b1978ae7752c41306e131a684e1f06"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "1.0.7"
+
+[[Pkg]]
+deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[SCS]]
+deps = ["BinaryProvider", "Libdl", "LinearAlgebra", "MathOptInterface", "MathProgBase", "SCS_jll", "SparseArrays"]
+git-tree-sha1 = "0b8dd347e857fe5cbb860df75169912e56e5ce2b"
+repo-rev = "master"
+repo-url = "https://github.com/jump-dev/SCS.jl.git"
+uuid = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
+version = "0.7.0"
+
+[[SCS_jll]]
+deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
+git-tree-sha1 = "293589b362649951e7dbe43f60f8b75ce0fae8fc"
+uuid = "f4f2fc5b-1d94-523c-97ea-2ab488bedf4b"
+version = "2.1.1+0"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[TimerOutputs]]
+deps = ["Printf"]
+git-tree-sha1 = "f458ca23ff80e46a630922c555d838303e4b9603"
+uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+version = "0.5.6"
+
+[[TranscodingStreams]]
+deps = ["Random", "Test"]
+git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.9.5"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[ZipFile]]
+deps = ["Libdl", "Printf", "Zlib_jll"]
+git-tree-sha1 = "254975fef2fc526583bb9b7c9420fe66ffe09f2f"
+uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
+version = "0.9.2"
+
+[[Zlib_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71"
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+version = "1.2.11+14"
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..e38e312
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
diff --git a/pmedian.jl b/pmedian.jl
new file mode 100644
index 0000000..177fb5b
--- /dev/null
+++ b/pmedian.jl
@@ -0,0 +1,271 @@
+import MathOptInterface
+import SCS
+import GLPK
+import Random
+
+using SparseArrays
+using TimerOutputs
+
+const MOI = MathOptInterface
+
+struct PMedianData
+    num_facilities::Int
+    num_customers::Int
+    num_locations::Int
+    customer_locations::Vector{Float64}
+end
+
+# This is the LP relaxation.
+function generate_moi_problem(model, data::PMedianData)
+    facility_variables = MOI.add_variables(model, data.num_locations)
+    for v in facility_variables
+        MOI.add_constraint(model, v, MOI.GreaterThan(0.0))
+        MOI.add_constraint(model, v, MOI.LessThan(1.0))
+    end
+    assignment_variables = [MOI.add_variable(model) for i in 1:data.num_customers, j in 1:data.num_locations]
+    for v in assignment_variables
+        MOI.add_constraint(model, v, MOI.GreaterThan(0.0))
+        # "Less than 1.0" constraint is redundant.
+    end
+    objective = MOI.ScalarAffineFunction(
+        [MOI.ScalarAffineTerm(abs(data.customer_locations[i] - j), assignment_variables[i,j])
+        for i in 1:data.num_customers for j in 1:data.num_locations],
+            0.0)
+    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    for i in 1:data.num_customers, j in 1:data.num_locations
+        # assignment_variables[i,j] <= facility_variables[j]
+        MOI.add_constraint(model, MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(1.0, assignment_variables[i,j]),
+            MOI.ScalarAffineTerm(-1.0, facility_variables[j])], 0.0),
+            MOI.LessThan(0.0))
+    end
+    for i in 1:data.num_customers
+        # sum_j assignment_variables[i,j] = 1
+        MOI.add_constraint(model, MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(1.0, assignment_variables[i,j]) for j in 1:data.num_locations],
+            0.0),
+            MOI.EqualTo(1.0)
+        )
+    end
+    # sum_j facility_variables[j] = num_facilities
+    MOI.add_constraint(model, MOI.ScalarAffineFunction(
+        MOI.ScalarAffineTerm.(1.0, facility_variables), 0.0),
+        MOI.EqualTo{Float64}(data.num_facilities)
+    )
+    return assignment_variables, facility_variables
+end
+
+function generate_glpk_problem(prob, data::PMedianData)
+    facility_variables = [GLPK.glp_add_cols(prob, 1) for i in 1:data.num_locations]
+    assignment_variables = [GLPK.glp_add_cols(prob, 1) for i in 1:data.num_customers, j in 1:data.num_locations]
+    for j in 1:data.num_locations
+        GLPK.glp_set_col_bnds(prob, facility_variables[j], GLPK.GLP_DB, 0.0, 1.0)
+    end
+    for i in 1:data.num_customers, j in 1:data.num_locations
+        GLPK.glp_set_col_bnds(prob, assignment_variables[i, j], GLPK.GLP_DB, 0.0, 1.0)
+        GLPK.glp_set_obj_coef(prob, assignment_variables[i, j], abs(data.customer_locations[i] - j))
+    end
+    GLPK.glp_set_obj_dir(prob, GLPK.GLP_MIN)
+
+    I = Cint[0] # Extra padding because GLPK starts reading from the second element
+    J = Cint[0]
+    V = Float64[0.0]
+    for i in 1:data.num_customers, j in 1:data.num_locations
+        # assignment_variables[i,j] <= facility_variables[j]
+        index = GLPK.glp_add_rows(prob, 1)
+        GLPK.glp_set_row_bnds(prob, index, GLPK.GLP_UP, 0.0, 0.0)
+        push!(I, index)
+        push!(J, assignment_variables[i, j])
+        push!(V, 1.0)
+        push!(I, index)
+        push!(J, facility_variables[j])
+        push!(V, -1.0)
+    end
+    for i in 1:data.num_customers
+        # sum_j assignment_variables[i,j] = 1
+        index = GLPK.glp_add_rows(prob, 1)
+        GLPK.glp_set_row_bnds(prob, index, GLPK.GLP_FX, 1.0, 1.0)
+        for j in 1:data.num_locations
+            push!(I, index)
+            push!(J, assignment_variables[i, j])
+            push!(V, 1.0)
+        end
+    end
+    # sum_j facility_variables[j] = num_facilities
+    index = GLPK.glp_add_rows(prob, 1)
+    GLPK.glp_set_row_bnds(prob, index, GLPK.GLP_FX, data.num_facilities, data.num_facilities)
+    for j in 1:data.num_locations
+        push!(I, index)
+        push!(J, facility_variables[j])
+        push!(V, 1.0)
+    end
+    GLPK.glp_load_matrix(prob, length(I) - 1, I, J, V)
+end
+
+function solve_glpk_direct(data::PMedianData; time_limit_sec=Inf)
+    @timeit "GLPK direct" begin
+        prob = GLPK.glp_create_prob()
+        param = GLPK.glp_smcp()
+        GLPK.glp_init_smcp(param)
+        param.msg_lev = GLPK.GLP_MSG_ERR
+        if isfinite(time_limit_sec)
+            param.tm_lim = 1000 * time_limit_sec
+        end
+
+        @timeit "generate" generate_glpk_problem(prob, data)
+        @timeit "solve" GLPK.glp_simplex(prob, param)
+        objval = GLPK.glp_get_obj_val(prob)
+        GLPK.glp_delete_prob(prob)
+    end
+    return objval
+end
+
+function generate_scs_problem(data::PMedianData)
+    # Equalities must be ordered first.
+    num_variables = 0
+    facility_variables = [num_variables += 1 for i in 1:data.num_locations]
+    assignment_variables = [num_variables += 1 for i in 1:data.num_customers, j in 1:data.num_locations]
+
+    b = Float64[]
+    I = Int[]
+    J = Int[]
+    V = Float64[]
+    
+    num_constraints = 0
+    for i in 1:data.num_customers
+        # sum_j assignment_variables[i,j] = 1
+        num_constraints += 1
+        push!(b, 1.0)
+        for j in 1:data.num_locations
+            push!(I, num_constraints)
+            push!(J, assignment_variables[i, j])
+            push!(V, 1.0)
+        end
+    end
+    # sum_j facility_variables[j] = num_facilities
+    num_constraints += 1
+    push!(b, data.num_facilities)
+    for j in 1:data.num_locations
+        push!(I, num_constraints)
+        push!(J, facility_variables[j])
+        push!(V, 1.0)
+    end
+
+    # Now inequality constraints in the form b - a'x >= 0.
+
+    for i in 1:data.num_customers, j in 1:data.num_locations
+        # assignment_variables[i,j] <= facility_variables[j]
+        num_constraints += 1
+        push!(b, 0.0)
+        push!(I, num_constraints)
+        push!(J, assignment_variables[i, j])
+        push!(V, 1.0)
+        push!(I, num_constraints)
+        push!(J, facility_variables[j])
+        push!(V, -1.0)
+    end
+
+    # All variables are nonnegative.
+    for v in 1:num_variables
+        num_constraints += 1
+        push!(b, 0.0)
+        push!(I, num_constraints)
+        push!(J, v)
+        push!(V, -1.0)
+    end
+
+    # facility_variables <= 1
+    for j in 1:data.num_locations
+        num_constraints += 1
+        push!(b, 1.0)
+        push!(I, num_constraints)
+        push!(J, facility_variables[j])
+        push!(V, 1.0)
+    end
+
+    c = zeros(num_variables)
+    for i in 1:data.num_customers, j in 1:data.num_locations
+        c[assignment_variables[i, j]] = abs(data.customer_locations[i] - j)
+    end
+    num_zero_cones = data.num_customers + 1
+    num_positive_orthants = num_constraints - num_zero_cones
+
+    A = sparse(I, J, V, num_constraints, num_variables)
+
+    return num_constraints, num_variables, A, b, c, num_zero_cones, num_positive_orthants,
+        Int[], Int[], 0, 0, Float64[]
+end
+
+function solve_scs_direct(data::PMedianData; max_iters)
+    @timeit "SCS direct" begin
+        @timeit "generate" scs_prob = generate_scs_problem(data)
+        @timeit "solve" solution = SCS.SCS_solve(SCS.IndirectSolver, 
+                                    scs_prob...; max_iters=max_iters,
+                                    acceleration_lookback=0, verbose=0)
+    end
+    objval = scs_prob[5]'*solution.x
+    return objval
+end
+
+function solve_moi(data::PMedianData, optimizer; params)
+    model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer(
+        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+        optimizer()), Float64)
+    for (param, value) in params
+        MOI.set(model, param, value)
+    end
+    @timeit "generate" x, y = generate_moi_problem(model, data)
+    @timeit "solve" MOI.optimize!(model)
+    return MOI.get(model, MOI.ObjectiveValue())
+end
+
+
+function solve_glpk_moi(data::PMedianData; time_limit_sec=Inf)
+    params = []
+    if isfinite(time_limit_sec)
+        push!(params, (MOI.TimeLimitSec(), time_limit_sec))
+    end
+    @timeit "GLPK MOI" solve_moi(data, GLPK.Optimizer, params=params)
+end
+
+function solve_scs_moi(data::PMedianData; max_iters)
+    params = [(MOI.RawParameter("max_iters"), max_iters),
+              (MOI.Silent(), true),
+              (MOI.RawParameter("acceleration_lookback"), 0)]
+    @timeit "SCS MOI" solve_moi(data, SCS.Optimizer, params=params)
+end
+
+
+
+function run_benchmark(;num_facilities, num_customers, num_locations,
+        time_limit_sec, max_iters)
+    Random.seed!(10)
+    reset_timer!()
+    data = PMedianData(num_facilities, num_customers, num_locations, rand(num_customers) .* num_locations)
+    
+    glpk_moi_obj = solve_glpk_moi(data, time_limit_sec=time_limit_sec)
+    @show glpk_moi_obj
+
+    glpk_direct_obj = solve_glpk_direct(data, time_limit_sec=time_limit_sec)
+    @show glpk_direct_obj
+
+    scs_moi_obj = solve_scs_moi(data, max_iters=max_iters)
+    @show scs_moi_obj
+
+    scs_direct_obj = solve_scs_direct(data, max_iters=max_iters)
+    @show scs_direct_obj
+    
+    print_timer()
+    println()
+end
+
+# JIT warm-up
+run_benchmark(num_facilities=5, num_customers=20, num_locations=10,
+    time_limit_sec=Inf, max_iters=10000)
+
+run_benchmark(num_facilities=5, num_customers=20, num_locations=10,
+    time_limit_sec=Inf, max_iters=10000)
+
+run_benchmark(num_facilities=10, num_customers=2000, num_locations=1000,
+    time_limit_sec=5, max_iters=10)