diff --git a/.gitignore b/.gitignore
index 06c7e4945..b3de0c16b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@
 .ipynb_checkpoints/
 Manifest.toml
 benchmark/*.json
+test/Project.toml
diff --git a/Project.toml b/Project.toml
index 56119ffbd..a82a3fb26 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,7 @@ version = "0.12.6"
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -14,13 +14,13 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 AbstractTrees = "^0.2.1"
 BenchmarkTools = "^0.4"
-MathProgBase = "^0.7"
+MathOptInterface = "~0.9"
 OrderedCollections = "^1.0"
 julia = "^1.0"
 
 [extras]
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
-GLPKMathProgInterface = "3c7084bd-78ad-589a-b5bb-dbd673274bea"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
@@ -28,4 +28,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ECOS", "GLPKMathProgInterface", "LinearAlgebra", "Random", "SCS", "Statistics", "Test"]
+test = ["ECOS", "GLPK", "LinearAlgebra", "Random", "SCS", "Statistics", "Test"]
diff --git a/README.md b/README.md
index 957936ac6..04b31538b 100644
--- a/README.md
+++ b/README.md
@@ -5,7 +5,7 @@
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://www.juliaopt.org/Convex.jl/stable)
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://www.juliaopt.org/Convex.jl/dev)
 
-**Convex.jl** is a [Julia](http://julialang.org) package for [Disciplined Convex Programming](http://dcp.stanford.edu/). Convex.jl can solve linear programs, mixed-integer linear programs, and DCP-compliant convex programs using a variety of solvers, including [Mosek](https://github.com/JuliaOpt/Mosek.jl), [Gurobi](https://github.com/JuliaOpt/Gurobi.jl), [ECOS](https://github.com/JuliaOpt/ECOS.jl), [SCS](https://github.com/JuliaOpt/SCS.jl), and  [GLPK](https://github.com/JuliaOpt/GLPK.jl), through the [MathProgBase](http://mathprogbasejl.readthedocs.org/en/latest/) interface. It also supports optimization with complex variables and coefficients.
+**Convex.jl** is a [Julia](http://julialang.org) package for [Disciplined Convex Programming](http://dcp.stanford.edu/). Convex.jl can solve linear programs, mixed-integer linear programs, and DCP-compliant convex programs using a variety of solvers, including [Mosek](https://github.com/JuliaOpt/Mosek.jl), [Gurobi](https://github.com/JuliaOpt/Gurobi.jl), [ECOS](https://github.com/JuliaOpt/ECOS.jl), [SCS](https://github.com/JuliaOpt/SCS.jl), and  [GLPK](https://github.com/JuliaOpt/GLPK.jl), through [MathOptInterface](https://github.com/JuliaOpt/MathOptInterface.jl). It also supports optimization with complex variables and coefficients.
 
 **Installation**: `julia> Pkg.add("Convex")`
 
@@ -38,7 +38,7 @@ x = Variable(n)
 problem = minimize(sumsquares(A * x - b), [x >= 0])
 
 # Solve the problem by calling solve!
-solve!(problem, SCSSolver())
+solve!(problem, SCS.Optimizer())
 
 # Check the status of the problem
 problem.status # :Optimal, :Infeasible, :Unbounded etc.
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index c56b4dc86..16a1b2d97 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -2,12 +2,14 @@ using Pkg
 tempdir = mktempdir()
 Pkg.activate(tempdir)
 Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))
-Pkg.add(["BenchmarkTools", "PkgBenchmark"])
+Pkg.add(["BenchmarkTools", "PkgBenchmark", "MathOptInterface"])
 Pkg.resolve()
 
 using Convex: Convex, ProblemDepot
 using BenchmarkTools
-
+using MathOptInterface
+const MOI = MathOptInterface
+const MOIU = MOI.Utilities
 
 const SUITE = BenchmarkGroup()
 
@@ -30,4 +32,7 @@ problems =  [
                 "mip_integer_variables",
             ]
 
-SUITE["formulation"] = ProblemDepot.benchmark_suite(Convex.conic_problem, problems)
+SUITE["formulation"] = ProblemDepot.benchmark_suite(problems) do problem
+    model = MOIU.MockOptimizer(MOIU.Model{Float64}())
+    Convex.load_MOI_model!(model, problem)
+end
diff --git a/docs/Project.toml b/docs/Project.toml
index 1712b9532..c6d03ec9d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,7 +6,7 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
-GLPKMathProgInterface = "3c7084bd-78ad-589a-b5bb-dbd673274bea"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 Interact = "c601a237-2ae4-5e1e-952c-7a85b0c7eef1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
diff --git a/docs/examples_literate/general_examples/basic_usage.jl b/docs/examples_literate/general_examples/basic_usage.jl
index da5c89191..cd0065ee4 100644
--- a/docs/examples_literate/general_examples/basic_usage.jl
+++ b/docs/examples_literate/general_examples/basic_usage.jl
@@ -8,7 +8,7 @@ end
 
 using SCS
 ## passing in verbose=0 to hide output from SCS
-solver = SCSSolver(verbose=0)
+solver = SCS.Optimizer(verbose=0)
 
 # ### Linear program
 #
@@ -73,7 +73,7 @@ p.optval
 
 x = Variable(4)
 p = satisfy(norm(x) <= 100, exp(x[1]) <= 5, x[2] >= 7, geomean(x[3], x[4]) >= x[2])
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 println(p.status)
 x.value
 
@@ -82,7 +82,7 @@ x.value
 
 y = Semidefinite(2)
 p = maximize(lambdamin(y), tr(y)<=6)
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 p.optval
 
 #-
@@ -105,10 +105,10 @@ y.value
 # $$
 #
 
-using GLPKMathProgInterface
+using GLPK
 x = Variable(4, :Int)
 p = minimize(sum(x), x >= 0.5)
-solve!(p, GLPKSolverMIP())
+solve!(p, GLPK.Optimizer())
 x.value
 
 #-
diff --git a/docs/examples_literate/general_examples/chebyshev_center.jl b/docs/examples_literate/general_examples/chebyshev_center.jl
index e1ce13bbf..9175e6205 100644
--- a/docs/examples_literate/general_examples/chebyshev_center.jl
+++ b/docs/examples_literate/general_examples/chebyshev_center.jl
@@ -25,7 +25,7 @@ p.constraints += a1' * x_c + r * norm(a1, 2) <= b[1];
 p.constraints += a2' * x_c + r * norm(a2, 2) <= b[2];
 p.constraints += a3' * x_c + r * norm(a3, 2) <= b[3];
 p.constraints += a4' * x_c + r * norm(a4, 2) <= b[4];
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 p.optval
 
 # Generate the figure
diff --git a/docs/examples_literate/general_examples/control.jl b/docs/examples_literate/general_examples/control.jl
index 7b381c8ad..65dac2cbe 100644
--- a/docs/examples_literate/general_examples/control.jl
+++ b/docs/examples_literate/general_examples/control.jl
@@ -121,7 +121,7 @@ push!(constraints, velocity[:, T] == 0)
 
 ## Solve the problem
 problem = minimize(sumsquares(force), constraints)
-solve!(problem, SCSSolver(verbose=0))
+solve!(problem, SCS.Optimizer(verbose=0))
 
 # We can plot the trajectory taken by the object.
 
diff --git a/docs/examples_literate/general_examples/huber_regression.jl b/docs/examples_literate/general_examples/huber_regression.jl
index a15adc2ad..b2f0c1e0f 100644
--- a/docs/examples_literate/general_examples/huber_regression.jl
+++ b/docs/examples_literate/general_examples/huber_regression.jl
@@ -40,18 +40,18 @@ for i=1:length(p_vals)
     fit = norm(beta - beta_true) / norm(beta_true);
     cost = norm(X' * beta - Y);
     prob = minimize(cost);
-    solve!(prob, SCSSolver(verbose=0));
+    solve!(prob, SCS.Optimizer(verbose=0));
     lsq_data[i] = evaluate(fit);
     
     ## Form and solve a prescient regression problem,
     ## i.e., where the sign changes are known.
     cost = norm(factor .* (X'*beta) - Y);
-    solve!(minimize(cost), SCSSolver(verbose=0))
+    solve!(minimize(cost), SCS.Optimizer(verbose=0))
     prescient_data[i] = evaluate(fit);
     
     ## Form and solve the Huber regression problem.
     cost = sum(huber(X' * beta - Y, 1));
-    solve!(minimize(cost), SCSSolver(verbose=0))
+    solve!(minimize(cost), SCS.Optimizer(verbose=0))
     huber_data[i] = evaluate(fit);
 end
 
diff --git a/docs/examples_literate/general_examples/logistic_regression.jl b/docs/examples_literate/general_examples/logistic_regression.jl
index d1caba41c..c63a726af 100644
--- a/docs/examples_literate/general_examples/logistic_regression.jl
+++ b/docs/examples_literate/general_examples/logistic_regression.jl
@@ -23,7 +23,7 @@ X = hcat(ones(size(iris, 1)), iris.SepalLength, iris.SepalWidth, iris.PetalLengt
 n, p = size(X)
 beta = Variable(p)
 problem = minimize(logisticloss(-Y.*(X*beta)))
-solve!(problem, SCSSolver(verbose=false))
+solve!(problem, SCS.Optimizer(verbose=false))
 
 # Let's see how well the model fits.
 using Plots
diff --git a/docs/examples_literate/general_examples/max_entropy.jl b/docs/examples_literate/general_examples/max_entropy.jl
index c7a5e3e21..6a9dfb216 100644
--- a/docs/examples_literate/general_examples/max_entropy.jl
+++ b/docs/examples_literate/general_examples/max_entropy.jl
@@ -23,7 +23,7 @@ b = rand(m, 1);
 
 x = Variable(n);
 problem = maximize(entropy(x), sum(x) == 1, A * x <= b)
-solve!(problem, SCSSolver(verbose=false))
+solve!(problem, SCS.Optimizer(verbose=false))
 problem.optval
 
 #-
diff --git a/docs/examples_literate/general_examples/optimal_advertising.jl b/docs/examples_literate/general_examples/optimal_advertising.jl
index a47af4f37..6a0997970 100644
--- a/docs/examples_literate/general_examples/optimal_advertising.jl
+++ b/docs/examples_literate/general_examples/optimal_advertising.jl
@@ -46,7 +46,7 @@ D = Variable(m, n);
 Si = [min(R[i]*dot(P[i,:], D[i,:]'), B[i]) for i=1:m];
 problem = maximize(sum(Si),
                [D >= 0, sum(D, dims=1)' <= T, sum(D, dims=2) >= c]);
-solve!(problem, SCSSolver(verbose=0));
+solve!(problem, SCS.Optimizer(verbose=0));
 
 #-
 
diff --git a/docs/examples_literate/general_examples/robust_approx_fitting.jl b/docs/examples_literate/general_examples/robust_approx_fitting.jl
index 01ca3fd64..7de9545a2 100644
--- a/docs/examples_literate/general_examples/robust_approx_fitting.jl
+++ b/docs/examples_literate/general_examples/robust_approx_fitting.jl
@@ -43,18 +43,18 @@ x = Variable(n)
 
 # Case 1: Nominal optimal solution
 p = minimize(norm(A * x - b, 2))
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 x_nom = evaluate(x)
 
 # Case 2: Stochastic robust approximation
 P = 1 / 3 * B' * B;
 p = minimize(square(pos(norm(A * x - b))) + quadform(x, Symmetric(P)))
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 x_stoch = evaluate(x)
 
 # Case 3: Worst-case robust approximation
 p = minimize(max(norm((A - B) * x - b), norm((A + B) * x - b)))
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 x_wc = evaluate(x)
 
 # Plot residuals:
diff --git a/docs/examples_literate/general_examples/svm.jl b/docs/examples_literate/general_examples/svm.jl
index 7d40068f8..c1662c9dd 100644
--- a/docs/examples_literate/general_examples/svm.jl
+++ b/docs/examples_literate/general_examples/svm.jl
@@ -34,7 +34,7 @@ neg_data = rand(MvNormal([-1.0, 2.0], 1.0), M);
 
 #-
 
-function svm(pos_data, neg_data, solver=SCSSolver(verbose=0))
+function svm(pos_data, neg_data, solver=SCS.Optimizer(verbose=0))
     ## Create variables for the separating hyperplane w'*x = b.
     w = Variable(n)
     b = Variable()
diff --git a/docs/examples_literate/general_examples/svm_l1regularization.jl b/docs/examples_literate/general_examples/svm_l1regularization.jl
index dcf22c77e..b64e8ffdc 100644
--- a/docs/examples_literate/general_examples/svm_l1regularization.jl
+++ b/docs/examples_literate/general_examples/svm_l1regularization.jl
@@ -34,8 +34,8 @@ beta_vals = zeros(length(beta), TRIALS);
 for i = 1:TRIALS
     lambda = lambda_vals[i];
     problem = minimize(loss/m + lambda*reg);
-    solve!(problem, ECOSSolver(verbose=0));
-    ## solve!(problem, SCSSolver(verbose=0,linear_solver=SCS.Direct, eps=1e-3))
+    solve!(problem, ECOS.Optimizer(verbose=0));
+    ## solve!(problem, SCS.Optimizer(verbose=0,linear_solver=SCS.Direct, eps=1e-3))
     train_error[i] = sum(float(sign.(X*beta_true .+ offset) .!= sign.(evaluate(X*beta - v))))/m;
     test_error[i] = sum(float(sign.(X_test*beta_true .+ offset) .!= sign.(evaluate(X_test*beta - v))))/TEST;
     beta_vals[:, i] =  evaluate(beta);
diff --git a/docs/examples_literate/general_examples/trade_off_curves.jl b/docs/examples_literate/general_examples/trade_off_curves.jl
index 17ff7bc1b..24f263176 100644
--- a/docs/examples_literate/general_examples/trade_off_curves.jl
+++ b/docs/examples_literate/general_examples/trade_off_curves.jl
@@ -18,7 +18,7 @@ x = Variable(n);
 for i=1:length(gammas)
     cost = sumsquares(A*x - b) + gammas[i]*norm(x,1);
     problem = minimize(cost, [norm(x, Inf) <= 1]);
-    solve!(problem, SCSSolver(verbose=0));
+    solve!(problem, SCS.Optimizer(verbose=0));
     x_values[:,i] = evaluate(x);
 end
 
diff --git a/docs/examples_literate/general_examples/worst_case_analysis.jl b/docs/examples_literate/general_examples/worst_case_analysis.jl
index 7725a0c01..4476fad5d 100644
--- a/docs/examples_literate/general_examples/worst_case_analysis.jl
+++ b/docs/examples_literate/general_examples/worst_case_analysis.jl
@@ -18,7 +18,7 @@ w = Variable(n);
 ret = dot(r, w);
 risk = sum(quadform(w, Sigma_nom));
 problem = minimize(risk, [sum(w) == 1, ret >= 0.1, norm(w, 1) <= 2])
-solve!(problem, SCSSolver(verbose=0));
+solve!(problem, SCS.Optimizer(verbose=0));
 wval = vec(evaluate(w))
 
 #-
@@ -31,7 +31,7 @@ problem = maximize(risk, [Sigma == Sigma_nom + Delta,
                     diag(Delta) == 0,
                     abs(Delta) <= 0.2,
                     Delta == Delta']);
-solve!(problem, SCSSolver(verbose=0));
+solve!(problem, SCS.Optimizer(verbose=0));
 println("standard deviation = ", round(sqrt(wval' * Sigma_nom * wval), sigdigits=2));
 println("worst-case standard deviation = ", round(sqrt(evaluate(risk)), sigdigits=2));
 println("worst-case Delta = ");
diff --git a/docs/examples_literate/mixed_integer/binary_knapsack.jl b/docs/examples_literate/mixed_integer/binary_knapsack.jl
index 10cd88137..97757186f 100644
--- a/docs/examples_literate/mixed_integer/binary_knapsack.jl
+++ b/docs/examples_literate/mixed_integer/binary_knapsack.jl
@@ -21,8 +21,8 @@ n = length(w)
 
 #-
 
-using Convex, GLPKMathProgInterface
+using Convex, GLPK
 x = Variable(n, :Bin)
 problem = maximize(dot(p, x), dot(w, x) <= C)
-solve!(problem, GLPKSolverMIP())
+solve!(problem, GLPK.Optimizer())
 evaluate(x)
diff --git a/docs/examples_literate/mixed_integer/n_queens.jl b/docs/examples_literate/mixed_integer/n_queens.jl
index 3f1d0cfbe..d2f45abd4 100644
--- a/docs/examples_literate/mixed_integer/n_queens.jl
+++ b/docs/examples_literate/mixed_integer/n_queens.jl
@@ -1,6 +1,6 @@
 # # N queens
 
-using Convex, GLPKMathProgInterface, LinearAlgebra, SparseArrays, Test
+using Convex, GLPK, LinearAlgebra, SparseArrays, Test
 aux(str) = joinpath(@__DIR__, "aux", str) # path to auxiliary files
 include(aux("antidiag.jl"))
 
@@ -16,7 +16,7 @@ constr += Constraint[sum(diag(x, k)) <= 1 for k = -n+2:n-2]
 ## Exactly one queen per row and one queen per column
 constr += Constraint[sum(x, dims=1) == 1, sum(x, dims=2) == 1]
 p = satisfy(constr)
-solve!(p, GLPKSolverMIP())
+solve!(p, GLPK.Optimizer())
 
 # Let us test the results:
 for k = -n+2:n-2
diff --git a/docs/examples_literate/mixed_integer/section_allocation.jl b/docs/examples_literate/mixed_integer/section_allocation.jl
index bcb20dc34..37c7017de 100644
--- a/docs/examples_literate/mixed_integer/section_allocation.jl
+++ b/docs/examples_literate/mixed_integer/section_allocation.jl
@@ -9,7 +9,7 @@
 # The goal will be to get an allocation matrix $X$, where $X_{ij} = 1$ if
 # student $i$ is assigned to section $j$ and $0$ otherwise. 
 
-using Convex, GLPKMathProgInterface
+using Convex, GLPK
 aux(str) = joinpath(@__DIR__, "aux", str) # path to auxiliary files
 
 # Load our preference matrix, `P`
@@ -30,5 +30,5 @@ constraints = [sum(X, dims=2) == 1, sum(X, dims=1) <= 10, sum(X, dims=1) >= 6]
 # students since the ranking of the first choice is 1.
 p = minimize(vec(X)' * vec(P), constraints)
 
-solve!(p, GLPKSolverMIP())
+solve!(p, GLPK.Optimizer())
 p.optval
diff --git a/docs/examples_literate/optimization_with_complex_variables/Fidelity in Quantum Information Theory.jl b/docs/examples_literate/optimization_with_complex_variables/Fidelity in Quantum Information Theory.jl
index 62f78dd96..e85aef53d 100644
--- a/docs/examples_literate/optimization_with_complex_variables/Fidelity in Quantum Information Theory.jl	
+++ b/docs/examples_literate/optimization_with_complex_variables/Fidelity in Quantum Information Theory.jl	
@@ -33,7 +33,7 @@ Z = ComplexVariable(n,n)
 objective = 0.5*real(tr(Z+Z'))
 constraint = [P Z;Z' Q] ⪰ 0
 problem = maximize(objective,constraint)
-solve!(problem, SCSSolver(verbose=0))
+solve!(problem, SCS.Optimizer(verbose=0))
 computed_fidelity = evaluate(objective)
 
 #-
diff --git a/docs/examples_literate/optimization_with_complex_variables/phase_recovery_using_MaxCut.jl b/docs/examples_literate/optimization_with_complex_variables/phase_recovery_using_MaxCut.jl
index 706c8a45c..d1bc404b9 100644
--- a/docs/examples_literate/optimization_with_complex_variables/phase_recovery_using_MaxCut.jl
+++ b/docs/examples_literate/optimization_with_complex_variables/phase_recovery_using_MaxCut.jl
@@ -69,7 +69,7 @@ objective = inner_product(U,M)
 c1 = diag(U) == 1 
 c2 = U in :SDP
 p = minimize(objective,c1,c2)
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 U.value
 
 #-
diff --git a/docs/examples_literate/optimization_with_complex_variables/povm_simulation.jl b/docs/examples_literate/optimization_with_complex_variables/povm_simulation.jl
index 7adbcd803..bc044d479 100644
--- a/docs/examples_literate/optimization_with_complex_variables/povm_simulation.jl
+++ b/docs/examples_literate/optimization_with_complex_variables/povm_simulation.jl
@@ -55,7 +55,7 @@ function get_visibility(K)
     constraints += t*K[3] + (1-t)*noise[3] == P[2][2] + P[4][2] + P[6][1]
     constraints += t*K[4] + (1-t)*noise[4] == P[3][2] + P[5][2] + P[6][2]
     p = maximize(t, constraints)
-    solve!(p, SCSSolver(verbose=0))
+    solve!(p, SCS.Optimizer(verbose=0))
     return p.optval
 end
 
diff --git a/docs/examples_literate/optimization_with_complex_variables/power_flow_optimization.jl b/docs/examples_literate/optimization_with_complex_variables/power_flow_optimization.jl
index da82e6b7f..06fe7bc9e 100644
--- a/docs/examples_literate/optimization_with_complex_variables/power_flow_optimization.jl
+++ b/docs/examples_literate/optimization_with_complex_variables/power_flow_optimization.jl
@@ -25,7 +25,7 @@ c3 = real(W[1,1]) == 1.06^2;
 push!(c1, c2)
 push!(c1, c3)
 p = maximize(objective, c1);
-solve!(p, SCSSolver(verbose = 0))
+solve!(p, SCS.Optimizer(verbose = 0))
 p.optval
 #15.125857662600703
 evaluate(objective)
diff --git a/docs/examples_literate/portfolio_optimization/portfolio_optimization.jl b/docs/examples_literate/portfolio_optimization/portfolio_optimization.jl
index 3fb5a90de..ca15785a3 100644
--- a/docs/examples_literate/portfolio_optimization/portfolio_optimization.jl
+++ b/docs/examples_literate/portfolio_optimization/portfolio_optimization.jl
@@ -56,7 +56,7 @@ p = minimize( risk,
               w_lower <= w, 
               w <= w_upper )
 
-solve!(p, SCSSolver())     #use SCSSolver(verbose = false) to suppress printing
+solve!(p, SCS.Optimizer())     #use SCS.Optimizer(verbose = false) to suppress printing
 
 #-
 
diff --git a/docs/examples_literate/portfolio_optimization/portfolio_optimization2.jl b/docs/examples_literate/portfolio_optimization/portfolio_optimization2.jl
index 177204771..37468be3c 100644
--- a/docs/examples_literate/portfolio_optimization/portfolio_optimization2.jl
+++ b/docs/examples_literate/portfolio_optimization/portfolio_optimization2.jl
@@ -52,7 +52,7 @@ for i = 1:N
     λ = λ_vals[i]
     p = minimize( λ*risk - (1-λ)*ret,
                   sum(w) == 1 )    
-    solve!(p, SCSSolver(verbose = false))
+    solve!(p, SCS.Optimizer(verbose = false))
     MeanVarA[i,:]= [evaluate(ret),evaluate(risk)[1]]    #risk is a 1x1 matrix
 end
 
@@ -68,7 +68,7 @@ for i = 1:N
                   sum(w) == 1,
                   w_lower <= w,     #w[i] is bounded
                   w <= w_upper )
-    solve!(p, SCSSolver(verbose = false))
+    solve!(p, SCS.Optimizer(verbose = false))
     MeanVarB[i,:]= [evaluate(ret),evaluate(risk)[1]]    
 end
 
@@ -95,7 +95,7 @@ for i = 1:N
     p = minimize( λ*risk - (1-λ)*ret,
                   sum(w) == 1,
                   (norm(w, 1)-1) <= Lmax)
-    solve!(p, SCSSolver(verbose = false))
+    solve!(p, SCS.Optimizer(verbose = false))
     MeanVarC[i,:]= [evaluate(ret),evaluate(risk)[1]]    
 end
 
diff --git a/docs/examples_literate/supplemental_material/Convex.jl_intro_ISMP2015.jl b/docs/examples_literate/supplemental_material/Convex.jl_intro_ISMP2015.jl
index 8a3cb002f..8a045282a 100644
--- a/docs/examples_literate/supplemental_material/Convex.jl_intro_ISMP2015.jl
+++ b/docs/examples_literate/supplemental_material/Convex.jl_intro_ISMP2015.jl
@@ -41,7 +41,7 @@ problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                    x >= 0)
 
 ## Solve the problem by calling solve!
-solve!(problem, SCSSolver(verbose=0))
+solve!(problem, SCS.Optimizer(verbose=0))
 
 println("problem status is ", problem.status) # :Optimal, :Infeasible, :Unbounded etc.
 println("optimal value is ", problem.optval)
@@ -54,7 +54,7 @@ using Interact, Plots
     global A
     problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                    x >= 0)
-    solve!(problem, SCSSolver(verbose=0))
+    solve!(problem, SCS.Optimizer(verbose=0))
     histogram(evaluate(x), xlims=(0,3.5), label="x")
 end
 
@@ -128,7 +128,7 @@ p = minimize(objective, constraint)
 #-
 
 ## solve the problem
-solve!(p, SCSSolver(verbose=0))
+solve!(p, SCS.Optimizer(verbose=0))
 p.status
 
 #-
@@ -152,7 +152,7 @@ evaluate(objective)
 # to solve problem using a different solver, just import the solver package and pass the solver to the `solve!` method: eg
 #
 #     using Mosek
-#     solve!(p, MosekSolver())
+#     solve!(p, Mosek.Optimizer())
 
 #-
 
@@ -173,9 +173,9 @@ x = Variable(n)
 μ = 1
 problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1), 
                    x >= 0)
-@time solve!(problem, SCSSolver(verbose=0))
+@time solve!(problem, SCS.Optimizer(verbose=0))
 λ = 1.5
-@time solve!(problem, SCSSolver(verbose=0), warmstart = true)
+@time solve!(problem, SCS.Optimizer(verbose=0), warmstart = true)
 
 # # DCP examples
 
diff --git a/docs/examples_literate/supplemental_material/paper_examples.jl b/docs/examples_literate/supplemental_material/paper_examples.jl
index 3a9f5ae16..6340605cb 100644
--- a/docs/examples_literate/supplemental_material/paper_examples.jl
+++ b/docs/examples_literate/supplemental_material/paper_examples.jl
@@ -13,7 +13,7 @@ e = 0;
   end
   p = minimize(e, x>=1);
 end
-@time solve!(p, ECOSSolver(verbose=0))
+@time solve!(p, ECOS.Optimizer(verbose=0))
 
 # Indexing.
 println("Indexing example")
@@ -26,7 +26,7 @@ e = 0;
   end
   p = minimize(e, x >= ones(1000, 1));
 end
-@time solve!(p, ECOSSolver(verbose=0))
+@time solve!(p, ECOS.Optimizer(verbose=0))
 
 # Matrix constraints.
 println("Matrix constraint example")
@@ -37,7 +37,7 @@ b = randn(p, n);
 @time begin
   p = minimize(norm(vec(X)), A * X == b);
 end
-@time solve!(p, ECOSSolver(verbose=0))
+@time solve!(p, ECOS.Optimizer(verbose=0))
 
 # Transpose.
 println("Transpose example")
@@ -46,12 +46,12 @@ A = randn(5, 5);
 @time begin
   p = minimize(norm2(X - A), X' == X);
 end
-@time solve!(p, ECOSSolver(verbose=0))
+@time solve!(p, ECOS.Optimizer(verbose=0))
 
 n = 3
 A = randn(n, n);
 #@time begin
   X = Variable(n, n);
   p = minimize(norm(vec(X' - A)), X[1,1] == 1);
-  solve!(p, ECOSSolver(verbose=0))
+  solve!(p, ECOS.Optimizer(verbose=0))
 #end
diff --git a/docs/examples_literate/time_series/time_series.jl b/docs/examples_literate/time_series/time_series.jl
index 85b834901..2ff931ec7 100644
--- a/docs/examples_literate/time_series/time_series.jl
+++ b/docs/examples_literate/time_series/time_series.jl
@@ -40,7 +40,7 @@ eq_constraints = [ yearly[i] == yearly[i - 365] for i in 365 + 1 : n ]
 smoothing = 100
 smooth_objective = sumsquares(yearly[1 : n - 1] - yearly[2 : n])
 problem = minimize(sumsquares(temps - yearly) + smoothing * smooth_objective, eq_constraints);
-solve!(problem, ECOSSolver(maxit=200, verbose=0))
+solve!(problem, ECOS.Optimizer(maxit=200, verbose=0))
 residuals = temps - evaluate(yearly)
 
 ## Plot smooth fit
@@ -82,7 +82,7 @@ end
 ## Solve autoregressive problem
 ar_coef = Variable(ar_len)
 problem = minimize(sumsquares(residuals_mat * ar_coef - residuals[ar_len + 1 : end]))
-solve!(problem, ECOSSolver(max_iters=200, verbose=0))
+solve!(problem, ECOS.Optimizer(max_iters=200, verbose=0))
 
 ## plot autoregressive fit of daily fluctuations for a few days
 ar_range = 1:145
diff --git a/docs/examples_literate/tomography/tomography.jl b/docs/examples_literate/tomography/tomography.jl
index 068ae0c9d..94c92d2df 100644
--- a/docs/examples_literate/tomography/tomography.jl
+++ b/docs/examples_literate/tomography/tomography.jl
@@ -62,7 +62,7 @@ pixel_colors = Variable(num_pixels)
 ## to reflect that, we minimize a norm
 objective = sumsquares(line_mat * pixel_colors - line_vals)
 problem = minimize(objective)
-solve!(problem, ECOSSolver(verbose=0))
+solve!(problem, ECOS.Optimizer(verbose=0))
 
 rows = zeros(img_size*img_size)
 cols = zeros(img_size*img_size)
diff --git a/docs/src/advanced.md b/docs/src/advanced.md
index 5473cbe2d..7e1267a97 100644
--- a/docs/src/advanced.md
+++ b/docs/src/advanced.md
@@ -13,7 +13,7 @@ using Convex, SCS
 x = Variable()
 constraint = x >= 0
 p = minimize(x, constraint)
-solve!(p, SCSSolver())
+solve!(p, SCS.Optimizer())
 
 # Get the dual value for the constraint
 p.constraints[1].dual
@@ -41,13 +41,13 @@ x = Variable(n)
 # first solve
 lambda = 100
 problem = minimize(sumsquares(y - x) + lambda * sumsquares(x - 10))
-@time solve!(problem, SCSSolver())
+@time solve!(problem, SCS.Optimizer())
 
 # now warmstart
 # if the solver takes advantage of warmstarts, 
 # this run will be faster
 lambda = 105
-@time solve!(problem, SCSSolver(), warmstart=true)
+@time solve!(problem, SCS.Optimizer(), warmstart=true)
 ```
 
 Fixing and freeing variables
@@ -89,12 +89,12 @@ for i=1:10
     # first solve for x
     # with y fixed, the problem is convex
     fix!(y)
-    solve!(problem, SCSSolver(), warmstart = i > 1 ? true : false)
+    solve!(problem, SCS.Optimizer(), warmstart = i > 1 ? true : false)
     free!(y)
 
     # now solve for y with x fixed at the previous solution
     fix!(x)
-    solve!(problem, SCSSolver(), warmstart = true)
+    solve!(problem, SCS.Optimizer(), warmstart = true)
     free!(x)
 end
 ```
diff --git a/docs/src/complex-domain_optimization.md b/docs/src/complex-domain_optimization.md
index 97c34541d..b290e33a6 100644
--- a/docs/src/complex-domain_optimization.md
+++ b/docs/src/complex-domain_optimization.md
@@ -106,7 +106,7 @@ constraints = [partialtrace(ρ, 1, [2; 2]) == [1 0; 0 0]
                tr(ρ) == 1
                ρ in :SDP]
 p = satisfy(constraints)
-solve!(p, SCSSolver(verbose=false))
+solve!(p, SCS.Optimizer(verbose=false))
 p.status
 ```
 
diff --git a/docs/src/contributing.md b/docs/src/contributing.md
index 73e46b168..95e429b43 100644
--- a/docs/src/contributing.md
+++ b/docs/src/contributing.md
@@ -73,15 +73,14 @@ Then read the conic form code:
  -   We define data structures for conic objectives and conic
      constraints, and simple ways of combining them, in
      [conic\_form.jl](https://github.com/JuliaOpt/Convex.jl/blob/master/src/conic_form.jl)
- -   We convert the internal conic form representation into the
-     [standard form for conic
-     solvers](http://mathprogbasejl.readthedocs.io/en/latest/conic.html)
-     in the function
-     [conic\_problem](https://github.com/JuliaOpt/Convex.jl/blob/master/src/problems.jl#L97).
+ -   We load the internal conic form representation into the
+     [MathOptInterface](https://github.com/JuliaOpt/MathOptInterface.jl)
+     model in the function
+     [load\_MOI\_model!](https://github.com/JuliaOpt/Convex.jl/blob/master/src/solution.jl#L151).
  -   We solve problems (that is, pass the standard form of the problem
      to a solver, and put the solution back into the values of the
      appropriate variables) in
-     [solution.jl](https://github.com/JuliaOpt/Convex.jl/blob/master/src/solution.jl#L8).
+     [solve!](https://github.com/JuliaOpt/Convex.jl/blob/master/src/solution.jl#L205).
 
 You're now armed and dangerous. Go ahead and open an issue (or comment
 on a previous one) if you can't figure something out, or submit a PR if
diff --git a/docs/src/credits.md b/docs/src/credits.md
index 8a821d007..ef4ac51dd 100644
--- a/docs/src/credits.md
+++ b/docs/src/credits.md
@@ -1,13 +1,17 @@
 Credits
 =======
 
-Convex.jl was developed and maintained by:
+Convex.jl was created, developed, and maintained by:
 
 -   [Jenny Hong](http://www.stanford.edu/~jyunhong/)
 -   [Karanveer Mohan](http://www.stanford.edu/~kvmohan/)
 -   [Madeleine Udell](http://www.stanford.edu/~udell/)
 -   [David Zeng](http://www.stanford.edu/~dzeng0/)
 
+Convex.jl is currently developed and maintained by the Julia
+community; see [Contributors](https://github.com/JuliaOpt/Convex.jl/graphs/contributors)
+for more.
+
 The Convex.jl developers also thank:
 
 -   the [JuliaOpt](http://www.juliaopt.org/) team: [Iain
diff --git a/docs/src/faq.md b/docs/src/faq.md
index 9e9c2ae67..80c133127 100644
--- a/docs/src/faq.md
+++ b/docs/src/faq.md
@@ -14,8 +14,8 @@ How does Convex.jl differ from JuMP?
 ------------------------------------
 
 Convex.jl and JuMP are both modelling languages for mathematical
-programming embedded in Julia, and both interface with solvers via the
-MathProgBase interface, so many of the same solvers are available in
+programming embedded in Julia, and both interface with solvers via
+MathOptInterface, so many of the same solvers are available in
 both. Convex.jl converts problems to a standard conic form. This
 approach requires (and certifies) that the problem is convex and DCP
 compliant, and guarantees global optimality of the resulting solution.
diff --git a/docs/src/index.md b/docs/src/index.md
index d29c9db27..a02c3a664 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -20,9 +20,8 @@ Convex.jl supports many solvers, including
 [Gurobi](https://github.com/JuliaOpt/gurobi.jl),
 [ECOS](https://github.com/JuliaOpt/ECOS.jl),
 [SCS](https://github.com/karanveerm/SCS.jl) and
-[GLPK](https://github.com/JuliaOpt/GLPK.jl), through the
-[MathProgBase](http://mathprogbasejl.readthedocs.org/en/latest/)
-interface.
+[GLPK](https://github.com/JuliaOpt/GLPK.jl), through
+[MathOptInterface](https://github.com/JuliaOpt/MathOptInterface.jl).
 
 Note that Convex.jl was previously called CVX.jl. This package is under
 active development; we welcome bug reports and feature requests. For
diff --git a/docs/src/installation.md b/docs/src/installation.md
index 90665de3d..f90acd9dd 100644
--- a/docs/src/installation.md
+++ b/docs/src/installation.md
@@ -18,6 +18,6 @@ Pkg.add("SCS")
 
 To solve certain problems such as mixed integer programming problems you
 will need to install another solver as well, such as
-[GLPK](https://github.com/JuliaOpt/GLPKMathProgInterface.jl). If you
+[GLPK](https://github.com/JuliaOpt/GLPK.jl). If you
 wish to use other solvers, please read the section on
 [Solvers](@ref).
diff --git a/docs/src/operations.md b/docs/src/operations.md
index 8f3decc85..d48810506 100644
--- a/docs/src/operations.md
+++ b/docs/src/operations.md
@@ -4,8 +4,7 @@ Operations
 Convex.jl currently supports the following functions. These functions
 may be composed according to the [DCP](http://dcp.stanford.edu)
 composition rules to form new convex, concave, or affine expressions.
-Convex.jl transforms each problem into an equivalent [cone
-program](http://mathprogbasejl.readthedocs.org/en/latest/conic.html) in
+Convex.jl transforms each problem into an equivalent conic program in
 order to pass the problem to a specialized solver. Depending on the
 types of functions used in the problem, the conic constraints may
 include linear, second-order, exponential, or semidefinite constraints,
diff --git a/docs/src/problem_depot.md b/docs/src/problem_depot.md
index 7b0f072cb..58e6c6244 100644
--- a/docs/src/problem_depot.md
+++ b/docs/src/problem_depot.md
@@ -11,7 +11,7 @@ For example, to test the solver SCS on all the problems of the depot except the
 using Convex, SCS, Test
 @testset "SCS" begin
     Convex.ProblemDepot.run_tests(; exclude=[r"mip"]) do p
-        solve!(p, SCSSolver(verbose=0, eps=1e-6))
+        solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
     end
 end
 ```
@@ -44,7 +44,7 @@ function affine_negate_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T})
 
 this should be the same for every problem, except for the name, which is a description of the problem. It should include what kind of atoms it uses (`affine` in this case), so that certain kinds of atoms can be ruled out by the `exclude` keyword to [`Convex.ProblemDepot.run_tests`](@ref) and [`Convex.ProblemDepot.benchmark_suite`](@ref); for example, many solvers cannot solve mixed-integer problems, so `mip` is included in the name of such problems.
 
-Then begins the body of the problem. It is setup like any other Convex.jl problem, only `handle_problem!` is called instead of `solve!`. This allows particular solvers to be used (via e.g. choosing `handle_problem! = p -> solve!(p, solver)`), or for any other function of the problem (e.g. `handle_problem! = p -> Convex.conic_problem(p)` which is used for benchmarking problem formulation speed.) Tests should be included and gated behind `if test` blocks, so that tests can be skipped for benchmarking, or in the case that the problem is not in fact solved during `handle_problem!`.
+Then begins the body of the problem. It is setup like any other Convex.jl problem, only `handle_problem!` is called instead of `solve!`. This allows particular solvers to be used (via e.g. choosing `handle_problem! = p -> solve!(p, solver)`), or for any other function of the problem. Tests should be included and gated behind `if test` blocks, so that tests can be skipped for benchmarking, or in the case that the problem is not in fact solved during `handle_problem!`.
 
 The fact that the problems may not be solved during `handle_problem!` brings with it a small complication: any command that assumes the problem has been solved should be behind an `if test` check. For example, in some of the problems, `real(x.value)` is used, for a variable `x`; perhaps as
 
diff --git a/docs/src/quick_tutorial.md b/docs/src/quick_tutorial.md
index d12be116f..67d1f5a3c 100644
--- a/docs/src/quick_tutorial.md
+++ b/docs/src/quick_tutorial.md
@@ -33,7 +33,7 @@ x = Variable(n)
 problem = minimize(sumsquares(A * x - b), [x >= 0])
 
 # Solve the problem by calling solve!
-solve!(problem, SCSSolver(verbose=false))
+solve!(problem, SCS.Optimizer(verbose=false))
 
 # Check the status of the problem
 problem.status # :Optimal, :Infeasible, :Unbounded etc.
diff --git a/docs/src/solvers.md b/docs/src/solvers.md
index 43ff0c764..be74d07d2 100644
--- a/docs/src/solvers.md
+++ b/docs/src/solvers.md
@@ -27,20 +27,19 @@ the documentation for that solver.
 To use a specific solver, you can use the following syntax
 
 ```julia
-solve!(p, GurobiSolver())
-solve!(p, MosekSolver())
-solve!(p, GLPKSolverMIP())
-solve!(p, GLPKSolverLP())
-solve!(p, ECOSSolver())
-solve!(p, SCSSolver())
+solve!(p, Gurobi.Optimizer())
+solve!(p, Mosek.Optimizer())
+solve!(p, GLPK.Optimizer())
+solve!(p, ECOS.Optimizer())
+solve!(p, SCS.Optimizer())
 ```
 
 (Of course, the solver must be installed first.) For example, we can use
 GLPK to solve a MILP
 
 ```julia
-using GLPKMathProgInterface
-solve!(p, GLPKSolverMIP())
+using GLPK
+solve!(p, GLPK.Optimizer())
 ```
 
 Many of the solvers also allow options to be passed in. More details can
@@ -51,7 +50,7 @@ in quiet mode), we can do so by
 
 ```julia
 using SCS
-solve!(p, SCSSolver(verbose=false))
+solve!(p, SCS.Optimizer(verbose=false))
 ```
 
 If we wish to increase the maximum number of iterations for ECOS or SCS,
@@ -59,9 +58,9 @@ we can do so by
 
 ```julia
 using ECOS
-solve!(p, ECOSSolver(maxit=10000))
+solve!(p, ECOS.Optimizer(maxit=10000))
 using SCS
-solve!(p, SCSSolver(max_iters=10000))
+solve!(p, SCS.Optimizer(max_iters=10000))
 ```
 
 To turn off the problem status warning issued by Convex when a solver is
@@ -70,5 +69,5 @@ not able to solve a problem to optimality, use the keyword argument
 solver parameters:
 
 ```julia
-solve!(p, SCSSolver(verbose=false), verbose=false)
+solve!(p, SCS.Optimizer(verbose=false), verbose=false)
 ```
diff --git a/docs/src/types.md b/docs/src/types.md
index 0a4c0a42c..1cde06b01 100644
--- a/docs/src/types.md
+++ b/docs/src/types.md
@@ -73,7 +73,7 @@ y = Variable()
 z = Variable()
 expr = x + y + z
 problem = minimize(expr, x >= 1, y >= x, 4 * z >= y)
-solve!(problem, SCSSolver())
+solve!(problem, SCS.Optimizer())
 
 # Once the problem is solved, we can call evaluate() on expr:
 evaluate(expr)
@@ -148,7 +148,7 @@ A problem can be solved by calling `solve!`
 solve!(problem, solver)
 ```
 
-passing a solver such as `SCSSolver()` from the package `SCS` as the
+passing a solver such as `SCS.Optimizer()` from the package `SCS` as the
 second argument. After the problem is solved, `problem.status` records
 the status returned by the optimization solver, and can be `:Optimal`,
 `:Infeasible`, `:Unbounded`, `:Indeterminate` or `:Error`. If the status
diff --git a/src/Convex.jl b/src/Convex.jl
index ffd0d899b..0245d4c0f 100644
--- a/src/Convex.jl
+++ b/src/Convex.jl
@@ -5,6 +5,11 @@ using OrderedCollections: OrderedDict
 using LinearAlgebra
 using SparseArrays
 
+using MathOptInterface
+const MOI = MathOptInterface
+const MOIU = MOI.Utilities
+const MOIB = MOI.Bridges
+
 global DEFAULT_SOLVER = nothing
 ### modeling framework
 include("dcp.jl")
@@ -20,7 +25,6 @@ include("constraints/signs_and_sets.jl")
 include("constraints/soc_constraints.jl")
 include("constraints/exp_constraints.jl")
 include("constraints/sdp_constraints.jl")
-include("solver_info.jl")
 include("problems.jl")
 include("solution.jl")
 
diff --git a/src/atoms/affine/conv.jl b/src/atoms/affine/conv.jl
index 8b02f651f..02ed3d835 100644
--- a/src/atoms/affine/conv.jl
+++ b/src/atoms/affine/conv.jl
@@ -11,7 +11,7 @@ function conv(x::Value, y::AbstractExpr)
     end
     m = length(x)
     n = y.size[1]
-    X = spzeros(m+n - 1, n)
+    X = spzeros(eltype(x), m + n - 1, n)
     for i = 1:n
         X[i:m+i-1, i] = x
     end
diff --git a/src/conic_form.jl b/src/conic_form.jl
index b5ba98dac..45bcc920b 100644
--- a/src/conic_form.jl
+++ b/src/conic_form.jl
@@ -90,7 +90,8 @@ end
 
 # A conic constraint is of the form [affine_expr1, affine_expr2, ..., affine_exprk] \in cone
 # we represent each affine expressions as a ConicObj
-# we represent the cone as a Symbol (defined in MathProgBase), like :SOC, :LP, etc
+# we represent the cone as a Symbol, like :SOC, :SDP, etc.
+# See the details of `get_MOI_set` for the full list of symbols used.
 # and we record the sizes of the affine expressions (XXX check...)
 # XXX might it be better to represent objs as a single ConicObj rather than an array of them?
 struct ConicConstr
diff --git a/src/constraints/sdp_constraints.jl b/src/constraints/sdp_constraints.jl
index 62218cb8a..3758f3990 100644
--- a/src/constraints/sdp_constraints.jl
+++ b/src/constraints/sdp_constraints.jl
@@ -46,24 +46,22 @@ function conic_form!(c::SDPConstraint, unique_conic_forms::UniqueConicForms)
         # the upper triangular part (not including diagonal)
         # and the corresponding entries in the lower triangular part, so
         # symmetry => c.child[upperpart]
-        # scale off-diagonal elements by sqrt(2)
-        rescale = sqrt(2)*tril(ones(n,n))
+        rescale = triu(ones(n,n))
         rescale[diagind(n, n)] .= 1.0
-        diagandlowerpart = findall(!iszero, vec(rescale))
+        diagandupperpart = findall(!iszero, vec(rescale))
         lowerpart = Vector{Int}(undef, div(n*(n-1),2))
         upperpart = Vector{Int}(undef, div(n*(n-1),2))
-        klower = 0
-        # diagandlowerpart in column-major order:
-        # ie the (1,1), (2,1), ..., (n,1), (2,2), (3,2), ...
-        # consider using  and find(triu(ones(3,3)))
+        kupper = 0
+        # diagandupperpart in column-major order:
+        # consider using find(triu(ones(3,3)))
         @inbounds for j = 1:n
             for i = j+1:n
-                klower += 1
-                upperpart[klower] = n*(i-1) + j # (j,i)th element
-                lowerpart[klower] = n*(j-1) + i # (i,j)th element
+                kupper += 1
+                upperpart[kupper] = n*(i-1) + j # (j,i)th element
+                lowerpart[kupper] =n*(j-1) + i # (i,j)th element
             end
         end
-        objective = conic_form!((broadcast(*,rescale,c.child))[diagandlowerpart], unique_conic_forms)
+        objective = conic_form!((broadcast(*,rescale,c.child))[diagandupperpart], unique_conic_forms)
         sdp_constraint = ConicConstr([objective], :SDP, [div(n*(n+1),2)])
         cache_conic_form!(unique_conic_forms, c, sdp_constraint)
         # make sure upper and lower triangular part match in the solution
diff --git a/src/problem_depot/problem_depot.jl b/src/problem_depot/problem_depot.jl
index 6bf0f89a5..2eca46e7e 100644
--- a/src/problem_depot/problem_depot.jl
+++ b/src/problem_depot/problem_depot.jl
@@ -3,7 +3,8 @@
 
 module ProblemDepot
 using BenchmarkTools, Test
-
+using MathOptInterface
+const MOI = MathOptInterface
 using Convex
 using LinearAlgebra
 using LinearAlgebra: eigen, I, opnorm
diff --git a/src/problem_depot/problems/affine.jl b/src/problem_depot/problems/affine.jl
index 4d197e026..5daa70096 100644
--- a/src/problem_depot/problems/affine.jl
+++ b/src/problem_depot/problems/affine.jl
@@ -1,6 +1,7 @@
 @add_problem affine function affine_negate_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable()
-    p = minimize(-x, [x <= 0])
+    p = minimize(-x, [x <= 0]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -22,7 +23,8 @@ end
 
 @add_problem affine function affine_multiply_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(1)
-    p = minimize(2.0 * x, [x >= 2, x <= 4])
+    p = minimize(2.0 * x, [x >= 2, x <= 4]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -34,7 +36,8 @@ end
 
     x = Variable(2)
     A = 1.5 * eye(2)
-    p = minimize([2 2] * x, [A * x >= [1.1; 1.1]])
+    p = minimize([2 2] * x, [A * x >= [1.1; 1.1]]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -51,7 +54,7 @@ end
     k = -y * [1.0, 2.0, 3.0]
     c = [y <= 3.0, y >= 0.0, x >= ones(3), k <= x, x <= z]
     o = 3 * y
-    p = Problem(:minimize, o, c)
+    p = Problem{T}(:minimize, o, c)
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -60,7 +63,7 @@ end
         @test p.optval ≈ 3 atol=atol rtol=rtol
     end
 
-    p = Problem(:minimize, o, c...)
+    p = Problem{T}(:minimize, o, c...)
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -71,7 +74,8 @@ end
 
     # Check #274
     x = ComplexVariable(2,2)
-    p = minimize( real( [1.0im, 0.0]' * x * [1.0im, 0.0] ), [ x == [1.0 0.0; 0.0 1.0] ])
+    p = minimize( real( [1.0im, 0.0]' * x * [1.0im, 0.0] ), [ x == [1.0 0.0; 0.0 1.0] ]; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test p.optval ≈ 1.0 atol=atol rtol=rtol
@@ -80,7 +84,8 @@ end
 
 @add_problem affine function affine_dot_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2)
-    p = minimize(dot([2.0; 2.0], x), x >= [1.1; 1.1])
+    p = minimize(dot([2.0; 2.0], x), x >= [1.1; 1.1]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -93,7 +98,8 @@ end
 
 @add_problem affine function affine_dot_atom_for_matrix_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2,2)
-    p = minimize(dot(fill(2.0, (2,2)), x), x >= 1.1)
+    p = minimize(dot(fill(2.0, (2,2)), x), x >= 1.1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -107,7 +113,8 @@ end
 @add_problem affine function affine_add_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(1)
     y = Variable(1)
-    p = minimize(x + y, [x >= 3, y >= 2])
+    p = minimize(x + y, [x >= 3, y >= 2]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -118,7 +125,8 @@ end
     end
 
     x = Variable(1)
-    p = minimize(x, [eye(2) + x >= eye(2)])
+    p = minimize(x, [eye(2) + x >= eye(2)]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -129,7 +137,8 @@ end
     end
 
     y = Variable()
-    p = minimize(y - 5, y >= -1)
+    p = minimize(y - 5, y >= -1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -143,7 +152,8 @@ end
 @add_problem affine function affine_transpose_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2)
     c = ones(2, 1)
-    p = minimize(x' * c, x >= 1)
+    p = minimize(x' * c, x >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -155,7 +165,8 @@ end
 
     X = Variable(2, 2)
     c = ones(2, 1)
-    p = minimize(c' * X' * c, [X >= ones(2, 2)])
+    p = minimize(c' * X' * c, [X >= ones(2, 2)]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -173,7 +184,7 @@ end
     c = ones(1, cols)
     d = ones(rows, 1)
     p = minimize(c * x' * d + d' * x * c' + (c * x''''' * d)',
-                [x' >= r_2, x >= r, x''' >= r_2, x'' >= r])
+                [x' >= r_2, x >= r, x''' >= r_2, x'' >= r]; numeric_type = T)
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -187,7 +198,8 @@ end
 
 @add_problem affine function affine_index_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2)
-    p = minimize(x[1] + x[2], [x >= 1])
+    p = minimize(x[1] + x[2], [x >= 1]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -199,7 +211,8 @@ end
 
     x = Variable(3)
     I = [true true false]
-    p = minimize(sum(x[I]), [x >= 1])
+    p = minimize(sum(x[I]), [x >= 1]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -215,7 +228,8 @@ end
     X = Variable(rows, cols)
     A = randn(rows, cols)
     c = rand(1, n)
-    p = minimize(c * X[1:n, 5:5+n-1]' * c', X >= A)
+    p = minimize(c * X[1:n, 5:5+n-1]' * c', X >= A; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -229,7 +243,8 @@ end
 
 @add_problem affine function affine_sum_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2,2)
-    p = minimize(sum(x), x>=1)
+    p = minimize(sum(x), x>=1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -240,7 +255,8 @@ end
     end
 
     x = Variable(2,2)
-    p = minimize(sum(x) - 2*x[1,1], x>=1, x[1,1]<=2)
+    p = minimize(sum(x) - 2*x[1,1], x>=1, x[1,1]<=2; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -252,7 +268,8 @@ end
 
     x = Variable(10)
     a = rand(10, 1)
-    p = maximize(sum(x[2:6]), x <= a)
+    p = maximize(sum(x[2:6]), x <= a; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -265,7 +282,8 @@ end
 
 @add_problem affine function affine_diag_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2,2)
-    p = minimize(sum(diag(x,1)), x >= 1)
+    p = minimize(sum(diag(x,1)), x >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -276,7 +294,8 @@ end
     end
 
     x = Variable(4, 4)
-    p = minimize(sum(diag(x)), x >= 2)
+    p = minimize(sum(diag(x)), x >= 2; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -289,7 +308,8 @@ end
 
 @add_problem affine function affine_trace_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2,2)
-    p = minimize(tr(x), x >= 1)
+    p = minimize(tr(x), x >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -302,7 +322,8 @@ end
 
 @add_problem affine function affine_dot_multiply_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
-    p = maximize(sum(dot(*)(x,[1,2,3])), x<=1)
+    p = maximize(sum(dot(*)(x,[1,2,3])), x<=1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -313,7 +334,8 @@ end
     end
 
     x = Variable(3, 3)
-    p = maximize(sum(dot(*)(x,eye(3))), x<=1)
+    p = maximize(sum(dot(*)(x,eye(3))), x<=1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -324,7 +346,8 @@ end
     end
 
     x = Variable(5, 5)
-    p = minimize(x[1, 1], dot(*)(3,x) >= 3)
+    p = minimize(x[1, 1], dot(*)(3,x) >= 3; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -335,7 +358,8 @@ end
     end
 
     x = Variable(3,1)
-    p = minimize(sum(dot(*)(ones(3,3), x)), x>=1)
+    p = minimize(sum(dot(*)(ones(3,3), x)), x>=1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -346,7 +370,8 @@ end
     end
 
     x = Variable(1,3)
-    p = minimize(sum(dot(*)(ones(3,3), x)), x>=1)
+    p = minimize(sum(dot(*)(ones(3,3), x)), x>=1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -357,7 +382,8 @@ end
     end
 
     x = Variable(1, 3, Positive())
-    p = maximize(sum(dot(/)(x,[1 2 3])), x<=1)
+    p = maximize(sum(dot(/)(x,[1 2 3])), x<=1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -379,7 +405,8 @@ end
     A = rand(2, 3)
     X = Variable(3, 2)
     c = rand()
-    p = minimize(sum(reshape(X, 2, 3) + A), X >= c)
+    p = minimize(sum(reshape(X, 2, 3) + A), X >= c; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -390,7 +417,8 @@ end
     end
 
     b = rand(6)
-    p = minimize(sum(vec(X) + b), X >= c)
+    p = minimize(sum(vec(X) + b), X >= c; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -404,7 +432,8 @@ end
     c = ones(16, 1)
     reshaped = reshape(x, 16, 1)
     a = collect(1:16)
-    p = minimize(c' * reshaped, reshaped >= a)
+    p = minimize(c' * reshaped, reshaped >= a; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -419,7 +448,8 @@ end
 @add_problem affine function affine_hcat_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(4, 4)
     y = Variable(4, 6)
-    p = maximize(sum(x) + sum([y fill(4.0, 4)]), [x y fill(2.0, (4, 2))] <= 2)
+    p = maximize(sum(x) + sum([y fill(4.0, 4)]), [x y fill(2.0, (4, 2))] <= 2; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -436,7 +466,8 @@ end
     y = Variable(4, 6)
 
     # TODO: fix dimension mismatch [y 4*eye(4); x -ones(4, 6)]
-    p = maximize(sum(x) + sum([y 4*eye(4); x -ones(4, 6)]), [x;y'] <= 2)
+    p = maximize(sum(x) + sum([y 4*eye(4); x -ones(4, 6)]), [x;y'] <= 2; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -456,7 +487,8 @@ end
     end
 
     x = Variable(4)
-    p = minimize(sum(Diagonal(x)), x == [1, 2, 3, 4])
+    p = minimize(sum(Diagonal(x)), x == [1, 2, 3, 4]; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -468,7 +500,8 @@ end
 
     x = Variable(3)
     c = [1; 2; 3]
-    p = minimize(c' * Diagonal(x) * c, x >= 1, sum(x) == 10)
+    p = minimize(c' * Diagonal(x) * c, x >= 1, sum(x) == 10; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -478,18 +511,19 @@ end
     end
 
     x = Variable(3)
-    p = minimize(sum(x), x >= 1, Diagonal(x)[1, 2] == 1)
-    output = handle_problem!(p)
+    p = minimize(sum(x), x >= 1, Diagonal(x)[1, 2] == 1; numeric_type = T)
+
+    handle_problem!(p)
     if test
-        @test output === nothing
-        @test p.status != :Optimal
+        @test p.status != MOI.OPTIMAL
     end
 end
 
 @add_problem affine function affine_conv_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
     h = [1, -1]
-    p = minimize(sum(conv(h, x)) + sum(x), x >= 1, x <= 2)
+    p = minimize(sum(conv(h, x)) + sum(x), x >= 1, x <= 2; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -501,7 +535,8 @@ end
 
     x = Variable(3)
     h = [1, -1]
-    p = minimize(sum(conv(x, h)) + sum(x), x >= 1, x <= 2)
+    p = minimize(sum(conv(x, h)) + sum(x), x >= 1, x <= 2; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -515,24 +550,27 @@ end
 
 @add_problem affine function affine_satisfy_problems(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable()
-    p = satisfy(x >= 0)
+    p = satisfy(x >= 0; numeric_type = T)
+
     add_constraints!(p, x >= 1)
     add_constraints!(p, [x >= -1, x <= 4])
     handle_problem!(p)
     if test
-        @test p.status == :Optimal
+        @test p.status == MOI.OPTIMAL
     end
 
-    p = satisfy([x >= 0, x >= 1, x <= 2])
+    p = satisfy([x >= 0, x >= 1, x <= 2]; numeric_type = T)
+
     handle_problem!(p)
     if test
-        @test p.status == :Optimal
+        @test p.status == MOI.OPTIMAL
     end
 
-    p = maximize(1, [x >= 1, x <= 2])
+    p = maximize(1, [x >= 1, x <= 2]; numeric_type = T)
+
     handle_problem!(p)
     if test
-        @test p.status == :Optimal
+        @test p.status == MOI.OPTIMAL
     end
 
     constr = x >= 0
@@ -540,55 +578,48 @@ end
     constr += x <= 10
     constr2 = x >= 0
     constr2 += [x >= 2, x <= 3] + constr
-    p = satisfy(constr)
+    p = satisfy(constr; numeric_type = T)
+
     handle_problem!(p)
     if test
-        @test p.status == :Optimal
+        @test p.status == MOI.OPTIMAL
     end
 end
 
-@add_problem affine function affine_dual(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+@add_problem affine function affine_dualvalue(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable()
-    p = minimize(x, x >= 0)
+    p = minimize(x, x >= 0; numeric_type = T)
+
     handle_problem!(p)
     if test
-        if p.solution.has_dual
-            println("Solution object has dual value, checking for dual correctness.")
-                @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol
-        end
+        @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol
     end
 
     x = Variable()
-    p = maximize(x, x <= 0)
+    p = maximize(x, x <= 0; numeric_type = T)
+
     handle_problem!(p)
     if test
-        if p.solution.has_dual
-        println("Solution object has dual value, checking for dual correctness.")
-            @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol
-        end
+        @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol
     end
 
     x = Variable()
-    p = minimize(x, x >= 0, x == 2)
+    p = minimize(x, x >= 0, x == 2; numeric_type = T)
+
     handle_problem!(p)
     if test
-        if p.solution.has_dual
-        println("Solution object has dual value, checking for dual correctness.")
-            @test p.constraints[1].dual ≈ 0 atol=atol rtol=rtol
-            @test abs.(p.constraints[2].dual) ≈ 1 atol=atol rtol=rtol
-        end
+        @test p.constraints[1].dual ≈ 0 atol=atol rtol=rtol
+        @test abs.(p.constraints[2].dual) ≈ 1 atol=atol rtol=rtol
     end
 
     x = Variable(2)
     A = 1.5 * eye(2)
-    p = minimize(dot([2.0; 2.0], x), [A * x >= [1.1; 1.1]])
+    p = minimize(dot([2.0; 2.0], x), [A * x >= [1.1; 1.1]]; numeric_type = T)
+
     handle_problem!(p)
     if test
-        if p.solution.has_dual
-        println("Solution object has dual value, checking for dual correctness.")
         dual = [4/3; 4/3]
-            @test all(abs.(p.constraints[1].dual - dual) .<= atol)
-        end
+        @test all(abs.(p.constraints[1].dual - dual) .<= atol)
     end
 end
 
@@ -607,9 +638,12 @@ end
     Rt2 = ComplexVariable(d,d)
     Rt3 = ComplexVariable(d,d)
     S = rand(ComplexF64,d,d)
-    handle_problem!(satisfy(partialtranspose(Rt1, 1, dims) == S ))
-    handle_problem!(satisfy(partialtranspose(Rt2, 2, dims) == S ))
-    handle_problem!(satisfy(partialtranspose(Rt3, 3, dims) == S ))
+    handle_problem!(satisfy(partialtranspose(Rt1, 1, dims) == S; numeric_type = T))
+
+    handle_problem!(satisfy(partialtranspose(Rt2, 2, dims) == S; numeric_type = T))
+
+    handle_problem!(satisfy(partialtranspose(Rt3, 3, dims) == S; numeric_type = T))
+
 
         
     if test
diff --git a/src/problem_depot/problems/benchmark.jl b/src/problem_depot/problems/benchmark.jl
index bcf2a295e..5371057f5 100644
--- a/src/problem_depot/problems/benchmark.jl
+++ b/src/problem_depot/problems/benchmark.jl
@@ -1,5 +1,6 @@
 @add_problem constraints_benchmark function LT_constraints(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = [Variable() for _ = 1:1000]
     for (i, xi) in enumerate(x)
         push!(p.constraints, xi <= 1.0 * i)
@@ -9,7 +10,8 @@
 end
 
 @add_problem constraints_benchmark function LT_constraint(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = Variable(1000)
     push!(p.constraints, x <= collect(1.0:1000.0))
     handle_problem!(p)
@@ -17,7 +19,8 @@ end
 end
 
 @add_problem constraints_benchmark function GT_constraints(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = [Variable() for _ = 1:1000]
     for (i, xi) in enumerate(x)
         push!(p.constraints, xi >= 1.0 * i)
@@ -27,7 +30,8 @@ end
 end
 
 @add_problem constraints_benchmark function GT_constraint(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = Variable(1000)
     push!(p.constraints, x >= collect(1.0:1000.0))
     handle_problem!(p)
@@ -36,7 +40,8 @@ end
 
 
 @add_problem constraints_benchmark function equality_constraints(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = [Variable() for _ = 1:1000]
     for (i, xi) in enumerate(x)
         push!(p.constraints, xi == 1.0 * i)
@@ -46,7 +51,8 @@ end
 end
 
 @add_problem constraints_benchmark function equality_constraint(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = Variable(1000)
     push!(p.constraints, x == collect(1.0:1000.0))
     handle_problem!(p)
@@ -54,7 +60,8 @@ end
 end
 
 @add_problem constraints_benchmark function sdp_constraint(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = Variable(44, 44) # 990 vectorized entries
     push!(p.constraints, x ⪰ 0)
     handle_problem!(p)
@@ -62,7 +69,8 @@ end
 end
 
 @add_problem constraints_benchmark function sdp_constraints(handle_problem!, args...)
-    p = satisfy()
+    p = satisfy(; numeric_type = T)
+
     x = [Variable(4, 4) for _ = 1:100] # 1000 total vectorized entries
     for v in x
         push!(p.constraints, v ⪰ 0)
diff --git a/src/problem_depot/problems/constant.jl b/src/problem_depot/problems/constant.jl
index 2f7bca287..cf66ad8b5 100644
--- a/src/problem_depot/problems/constant.jl
+++ b/src/problem_depot/problems/constant.jl
@@ -9,7 +9,8 @@
     β = Variable(5)
     β.value = ones(5)
 
-    problem = minimize(sum(c * β), [β >= 0])
+    problem = minimize(sum(c * β), [β >= 0]; numeric_type = T)
+
     handle_problem!(problem)
     if test
         @test problem.optval ≈ evaluate(sum(c * β)) atol=atol rtol=rtol
@@ -22,13 +23,15 @@ end
     x = Variable(2)
     y = Variable(2)
     fix!(x, [1 1]')
-    prob = minimize(y'*(x+[2 2]'), [y>=0])
+    prob = minimize(y'*(x+[2 2]'), [y>=0]; numeric_type = T)
+
     handle_problem!(prob)
     if test
         @test prob.optval ≈ 0.0 atol=atol rtol=rtol
     end
 
-    prob = minimize(x'*y, [y>=0])
+    prob = minimize(x'*y, [y>=0]; numeric_type = T)
+
     handle_problem!(prob)
     if test
         @test prob.optval ≈ 0.0 atol=atol rtol=rtol
@@ -39,7 +42,8 @@ end
     x = Variable()
     y = Variable()
     fix!(x, 1.0)
-    prob = minimize(y*x, [y >= x, x >= 0.5])
+    prob = minimize(y*x, [y >= x, x >= 0.5]; numeric_type = T)
+
     handle_problem!(prob)
     if test
         @test prob.optval ≈ 1.0 atol=atol rtol=rtol
@@ -63,7 +67,8 @@ end
     p = Variable()
     fix!(p, 1.0)
     x = Variable(2,2)
-    prob = minimize( tr(p*x), [ x >= 1 ])
+    prob = minimize( tr(p*x), [ x >= 1 ]; numeric_type = T)
+
     handle_problem!(prob)
     if test
         @test prob.optval ≈ 2.0 atol=atol rtol=rtol
@@ -75,7 +80,8 @@ end
     x = ComplexVariable()
     fix!(x, 1.0 + im*1.0)
     y = Variable()
-    prob = minimize( real(x*y), [ y >= .5, real(x) >= .5, imag(x) >= 0])
+    prob = minimize( real(x*y), [ y >= .5, real(x) >= .5, imag(x) >= 0]; numeric_type = T)
+
     handle_problem!(prob)
     if test
         @test prob.optval ≈ .5 atol=atol rtol=rtol
@@ -107,7 +113,8 @@ end
     x = ComplexVariable(5)
     fix!(x, ones(5) + im*ones(5))
     y = Variable()
-    prob = minimize( real(y*sum(x)), [ y >= .5, real(x) >= .5, imag(x) >= 0])
+    prob = minimize( real(y*sum(x)), [ y >= .5, real(x) >= .5, imag(x) >= 0]; numeric_type = T)
+
     handle_problem!(prob)
     if test
         @test prob.optval ≈ 2.5 atol=atol rtol=rtol
diff --git a/src/problem_depot/problems/exp.jl b/src/problem_depot/problems/exp.jl
index b74f85896..73b891ef0 100644
--- a/src/problem_depot/problems/exp.jl
+++ b/src/problem_depot/problems/exp.jl
@@ -1,6 +1,7 @@
 @add_problem exp function exp_exp_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable()
-    p = minimize(exp(y), y>=0)
+    p = minimize(exp(y), y>=0; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -11,7 +12,8 @@
     end
 
     y = Variable()
-    p = minimize(exp(y), y>=1)
+    p = minimize(exp(y), y>=1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -22,7 +24,8 @@
     end
 
     y = Variable(5)
-    p = minimize(sum(exp(y)), y>=0)
+    p = minimize(sum(exp(y)), y>=0; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -33,7 +36,8 @@
     end
 
     y = Variable(5)
-    p = minimize(sum(exp(y)), y>=0)
+    p = minimize(sum(exp(y)), y>=0; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -45,7 +49,8 @@ end
 
 @add_problem exp function exp_log_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable()
-    p = maximize(log(y), y<=1)
+    p = maximize(log(y), y<=1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -55,7 +60,8 @@ end
     end
 
     y = Variable()
-    p = maximize(log(y), y<=2)
+    p = maximize(log(y), y<=2; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -65,7 +71,8 @@ end
     end
 
     y = Variable()
-    p = maximize(log(y), [y<=2, exp(y)<=10])
+    p = maximize(log(y), [y<=2, exp(y)<=10]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -77,7 +84,8 @@ end
 
 @add_problem exp function exp_log_sum_exp_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable(5)
-    p = minimize(logsumexp(y), y>=1)
+    p = minimize(logsumexp(y), y>=1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -89,7 +97,8 @@ end
 
 @add_problem exp function exp_logistic_loss_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable(5)
-    p = minimize(logisticloss(y), y>=1)
+    p = minimize(logisticloss(y), y>=1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -101,7 +110,8 @@ end
 
 @add_problem exp function exp_entropy_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable(5, Positive())
-    p = maximize(entropy(y), sum(y)<=1)
+    p = maximize(entropy(y), sum(y)<=1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -115,7 +125,8 @@ end
     x = Variable(1)
     y = Variable(1)
     # x log (x/y)
-    p = minimize(relative_entropy(x,y), y==1, x >= 2)
+    p = minimize(relative_entropy(x,y), y==1, x >= 2; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -129,7 +140,8 @@ end
     x = Variable(1)
     y = Variable(1)
     # y log (x/y)
-    p = maximize(log_perspective(x,y), y==5, x <= 10)
+    p = maximize(log_perspective(x,y), y==5, x <= 10; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
diff --git a/src/problem_depot/problems/lp.jl b/src/problem_depot/problems/lp.jl
index 367039e98..caa407a9c 100644
--- a/src/problem_depot/problems/lp.jl
+++ b/src/problem_depot/problems/lp.jl
@@ -1,6 +1,7 @@
 @add_problem lp function lp_dual_abs_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable()
-    p = minimize(abs(x), x<=-1)
+    p = minimize(abs(x), x<=-1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -8,10 +9,12 @@
     if test
         @test p.optval ≈ 1 atol=atol rtol=rtol
         @test evaluate(abs(x)) ≈ 1 atol=atol rtol=rtol
+        @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol
     end
 
     x = Variable(2,2)
-    p = minimize(sum(abs(x)), x[2,2]>=1, x[1,1]>=1, x>=0)
+    p = minimize(sum(abs(x)), x[2,2]>=1, x[1,1]>=1, x>=0; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -30,7 +33,8 @@ end
 @add_problem lp function lp_maximum_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(10)
     a = shuffle(collect(0.1:0.1:1.0))
-    p = minimize(maximum(x), x >= a)
+    p = minimize(maximum(x), x >= a; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -44,7 +48,8 @@ end
 @add_problem lp function lp_minimum_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(1)
     a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
-    p = maximize(minimum(x), x <= a)
+    p = maximize(minimum(x), x <= a; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -61,7 +66,8 @@ end
     d = fill(2.0, (6, 1))
     constraints = [[x y] <= 2, z <= 0, z <= x, 2z >= -1]
     objective = sum(x + z) + minimum(y) + c' * y * d
-    p = maximize(objective, constraints)
+    p = maximize(objective, constraints; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -77,7 +83,8 @@ end
     y = Variable(10, 10)
     a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
     b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
-    p = minimize(maximum(max(x, y)), [x >= a, y >= b])
+    p = minimize(maximum(max(x, y)), [x >= a, y >= b]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -95,7 +102,8 @@ end
     y = Variable(10, 10)
     a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
     b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
-    p = maximize(minimum(min(x, y)), [x <= a, y <= b])
+    p = maximize(minimum(min(x, y)), [x <= a, y <= b]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -111,7 +119,8 @@ end
 @add_problem lp function lp_pos_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
     a = [-2; 1; 2]
-    p = minimize(sum(pos(x)), [x >= a, x <= 2])
+    p = minimize(sum(pos(x)), [x >= a, x <= 2]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -124,7 +133,8 @@ end
 
 @add_problem lp function lp_neg_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
-    p = minimize(1, [x >= -2, x <= -2, neg(x) <= 3])
+    p = minimize(1, [x >= -2, x <= -2, neg(x) <= 3]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -137,7 +147,8 @@ end
 
 @add_problem lp function lp_sumlargest_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2)
-    p = minimize(sumlargest(x, 2), x >= [1; 1])
+    p = minimize(sumlargest(x, 2), x >= [1; 1]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -148,7 +159,8 @@ end
     end
 
     x = Variable(4, 4)
-    p = minimize(sumlargest(x, 3), x >= eye(4), x[1, 1] >= 1.5, x[2, 3] >= 2.1)
+    p = minimize(sumlargest(x, 3), x >= eye(4), x[1, 1] >= 1.5, x[2, 3] >= 2.1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -161,7 +173,8 @@ end
 
 @add_problem lp function lp_sumsmallest_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(4, 4)
-    p = minimize(sumlargest(x, 2), sumsmallest(x, 4) >= 1)
+    p = minimize(sumlargest(x, 2), sumsmallest(x, 4) >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -172,7 +185,8 @@ end
     end
 
     x = Variable(3, 2)
-    p = maximize(sumsmallest(x, 3), x >= 2, x <= 5, sumlargest(x, 3) <= 12)
+    p = maximize(sumsmallest(x, 3), x >= 2, x <= 5, sumlargest(x, 3) <= 12; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -185,7 +199,8 @@ end
 
 @add_problem lp function lp_dotsort_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(4, 1)
-    p = minimize(dotsort(x, [1, 2, 3, 4]), sum(x) >= 7, x >= 0, x <= 2, x[4] <= 1)
+    p = minimize(dotsort(x, [1, 2, 3, 4]), sum(x) >= 7, x >= 0, x <= 2, x[4] <= 1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -197,7 +212,8 @@ end
     end
 
     x = Variable(2, 2)
-    p = minimize(dotsort(x, [1 2; 3 4]), sum(x) >= 7, x >= 0, x <= 2, x[2, 2] <= 1)
+    p = minimize(dotsort(x, [1 2; 3 4]), sum(x) >= 7, x >= 0, x <= 2, x[2, 2] <= 1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -214,7 +230,8 @@ end
 
 @add_problem lp function lp_dual_norm_inf_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
-    p = minimize(norm_inf(x), [-2 <= x, x <= 1])
+    p = minimize(norm_inf(x), [-2 <= x, x <= 1]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -229,7 +246,8 @@ end
 
 @add_problem lp function lp_dual_norm_1_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
-    p = minimize(norm_1(x), [-2 <= x, x <= 1])
+    p = minimize(norm_1(x), [-2 <= x, x <= 1]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
diff --git a/src/problem_depot/problems/mip.jl b/src/problem_depot/problems/mip.jl
index 2d454aaeb..118e952bd 100644
--- a/src/problem_depot/problems/mip.jl
+++ b/src/problem_depot/problems/mip.jl
@@ -1,6 +1,7 @@
 @add_problem mip function mip_lp_fallback_interface(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable()
-    p = minimize(x, x>=4.3)
+    p = minimize(x, x>=4.3; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -10,7 +11,8 @@
     end
 
     x = Variable(2)
-    p = minimize(norm(x,1), x[1]>=4.3)
+    p = minimize(norm(x,1), x[1]>=4.3; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -22,7 +24,8 @@ end
 
 @add_problem mip function mip_integer_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(:Int)
-    p = minimize(x, x>=4.3)
+    p = minimize(x, x>=4.3; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -32,7 +35,8 @@ end
     end
 
     x = Variable(2, :Int)
-    p = minimize(sum(x), x>=4.3)
+    p = minimize(sum(x), x>=4.3; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -43,7 +47,8 @@ end
 
     x = Variable(:Int)
     y = Variable()
-    p = minimize(sum(x + y), x>=4.3, y>=7)
+    p = minimize(sum(x + y), x>=4.3, y>=7; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -53,7 +58,8 @@ end
     end
 
     x = Variable(2, :Int)
-    p = minimize(norm(x, 1), x[1]>=4.3)
+    p = minimize(norm(x, 1), x[1]>=4.3; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -63,7 +69,8 @@ end
     end
 
     x = Variable(2, :Int)
-    p = minimize(sum(x), x[1]>=4.3, x>=0)
+    p = minimize(sum(x), x[1]>=4.3, x>=0; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -73,7 +80,8 @@ end
     end
 
     x = Variable(2, :Int)
-    p = minimize(sum(x), x>=.5)
+    p = minimize(sum(x), x>=.5; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -85,7 +93,8 @@ end
 
 @add_problem mip function mip_binary_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2, :Bin)
-    p = minimize(sum(x), x>=.5)
+    p = minimize(sum(x), x>=.5; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -95,7 +104,8 @@ end
     end
 
     x = Variable(2, :Bin)
-    p = minimize(sum(x), x[1]>=.5, x>=0)
+    p = minimize(sum(x), x[1]>=.5, x>=0; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
diff --git a/src/problem_depot/problems/sdp.jl b/src/problem_depot/problems/sdp.jl
index 0c84be4da..ea013f166 100644
--- a/src/problem_depot/problems/sdp.jl
+++ b/src/problem_depot/problems/sdp.jl
@@ -1,7 +1,8 @@
 # TODO: uncomment vexity checks once SDP on vars/constraints changes vexity of problem
 @add_problem sdp function sdp_sdp_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable((2,2), :Semidefinite)
-    p = minimize(y[1,1])
+    p = minimize(y[1,1]; numeric_type = T)
+
     # @fact vexity(p) --> ConvexVexity()
     handle_problem!(p)
     if test
@@ -9,7 +10,8 @@
     end
 
     y = Variable((3,3), :Semidefinite)
-    p = minimize(y[1,1], y[2,2]==1)
+    p = minimize(y[1,1], y[2,2]==1; numeric_type = T)
+
     # @fact vexity(p) --> ConvexVexity()
     handle_problem!(p)
     if test
@@ -20,13 +22,15 @@
     # This test fails on Mosek. See
     # https://github.com/JuliaOpt/Mosek.jl/issues/29
     # y = Variable((2, 2), :Semidefinite)
-    # p = minimize(y[1, 1], y[1, 2] == 1)
+    # p = minimize(y[1, 1], y[1, 2] == 1; numeric_type = T)
+
     # # @fact vexity(p) --> ConvexVexity()
     # handle_problem!(p)
     # @fact p.optval --> roughly(0, atol)
 
     y = Semidefinite(3)
-    p = minimize(sum(diag(y)), y[1, 1] == 1)
+    p = minimize(sum(diag(y)), y[1, 1] == 1; numeric_type = T)
+
     # @fact vexity(p) --> ConvexVexity()
     handle_problem!(p)
     if test
@@ -34,7 +38,8 @@
     end
 
     y = Variable((3, 3), :Semidefinite)
-    p = minimize(tr(y), y[2,1]<=4, y[2,2]>=3)
+    p = minimize(tr(y), y[2,1]<=4, y[2,2]>=3; numeric_type = T)
+
     # @fact vexity(p) --> ConvexVexity()
     handle_problem!(p)
     if test
@@ -43,7 +48,8 @@
 
     x = Variable(Positive())
     y = Semidefinite(3)
-    p = minimize(y[1, 2], y[2, 1] == 1)
+    p = minimize(y[1, 2], y[2, 1] == 1; numeric_type = T)
+
     # @fact vexity(p) --> ConvexVexity()
     handle_problem!(p)
     if test
@@ -55,7 +61,8 @@ end
     # This test fails on Mosek
     x = Variable(Positive())
     y = Variable((3, 3))
-    p = minimize(x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1)
+    p = minimize(x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1; numeric_type = T)
+
     # @fact vexity(p) --> ConvexVexity()
     handle_problem!(p)
     if test
@@ -65,7 +72,8 @@ end
 
 @add_problem sdp function sdp_nuclear_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Semidefinite(3)
-    p = minimize(nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2)
+    p = minimize(nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -78,7 +86,8 @@ end
 
 @add_problem sdp function sdp_operator_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable((3,3))
-    p = minimize(opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12)
+    p = minimize(opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -91,7 +100,8 @@ end
 
 @add_problem sdp function sdp_sigma_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable((3,3))
-    p = minimize(sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12)
+    p = minimize(sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -102,9 +112,10 @@ end
     end
 end
 
-@add_problem sdp function sdp_lambda_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+@add_problem sdp function sdp_dual_lambda_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Semidefinite(3)
-    p = minimize(lambdamax(y), y[1,1]>=4)
+    p = minimize(lambdamax(y), y[1,1]>=4; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -113,11 +124,22 @@ end
         @test p.optval ≈ 4 atol=atol rtol=rtol
         @test evaluate(lambdamax(y)) ≈ 4 atol=atol rtol=rtol
     end
+
+    # https://github.com/JuliaOpt/Convex.jl/issues/337
+    x = ComplexVariable(2, 2)
+    p = minimize( lambdamax(x), [ x[1,2] == im, x[2,2] == 1.0, x ⪰ - eye(2) ]; numeric_type = T)
+    handle_problem!(p)
+    if test
+        @test p.optval ≈ 1.5 atol=atol rtol=rtol
+        @test p.constraints[1].dual ≈ im atol=atol rtol=rtol
+        @test p.constraints[2].dual ≈ 0.75 atol=atol rtol=rtol
+    end
 end
 
 @add_problem sdp function sdp_lambda_min_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Semidefinite(3)
-    p = maximize(lambdamin(y), tr(y)<=6)
+    p = maximize(lambdamin(y), tr(y)<=6; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -131,7 +153,8 @@ end
 @add_problem sdp function sdp_matrix_frac_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = [1, 2, 3]
     P = Variable(3, 3)
-    p = minimize(matrixfrac(x, P), P <= 2*eye(3), P >= 0.5 * eye(3))
+    p = minimize(matrixfrac(x, P), P <= 2*eye(3), P >= 0.5 * eye(3); numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -145,7 +168,8 @@ end
 @add_problem sdp function sdp_matrix_frac_atom_both_arguments_variable(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
     P = Variable(3, 3)
-    p = minimize(matrixfrac(x, P), lambdamax(P) <= 2, x[1] >= 1)
+    p = minimize(matrixfrac(x, P), lambdamax(P) <= 2, x[1] >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -158,7 +182,8 @@ end
 
 @add_problem sdp function sdp_sum_largest_eigs(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Semidefinite(3)
-    p = minimize(sumlargesteigs(x, 2), x >= 1)
+    p = minimize(sumlargesteigs(x, 2), x >= 1; numeric_type = T)
+
     handle_problem!(p)
 
     if test
@@ -167,7 +192,8 @@ end
     end
 
     x = Semidefinite(3)
-    p = minimize(sumlargesteigs(x, 2), [x[i,:] >= i for i=1:3]...)
+    p = minimize(sumlargesteigs(x, 2), [x[i,:] >= i for i=1:3]...; numeric_type = T)
+
     handle_problem!(p)
 
     if test
@@ -175,11 +201,13 @@ end
     end
 
     x1 = Semidefinite(3)
-    p1 = minimize(lambdamax(x1), x1[1,1]>=4)
+    p1 = minimize(lambdamax(x1), x1[1,1]>=4; numeric_type = T)
+
     handle_problem!(p1)
 
     x2 = Semidefinite(3)
-    p2 = minimize(sumlargesteigs(x2, 1), x2[1,1]>=4)
+    p2 = minimize(sumlargesteigs(x2, 1), x2[1,1]>=4; numeric_type = T)
+
     handle_problem!(p2)
 
     if test
@@ -187,11 +215,13 @@ end
     end
 
     x1 = Semidefinite(3)
-    p1 = minimize(lambdamax(x1), [x1[i,:] >= i for i=1:3]...)
+    p1 = minimize(lambdamax(x1), [x1[i,:] >= i for i=1:3]...; numeric_type = T)
+
     handle_problem!(p1)
 
     x2 = Semidefinite(3)
-    p2 = minimize(sumlargesteigs(x2, 1), [x2[i,:] >= i for i=1:3]...)
+    p2 = minimize(sumlargesteigs(x2, 1), [x2[i,:] >= i for i=1:3]...; numeric_type = T)
+
     handle_problem!(p2)
 
     if test
@@ -204,7 +234,8 @@ end
     id = eye(4)
     X = Semidefinite(4)
     W = kron(id, X)
-    p = maximize(tr(W), tr(X) ≤ 1)
+    p = maximize(tr(W), tr(X) ≤ 1; numeric_type = T)
+
     if test
         @test vexity(p) == AffineVexity()
     end
@@ -219,7 +250,8 @@ end
     B = [1 0; 0 0]
     ρ = kron(B, A)
     constraints = [partialtrace(ρ, 1, [2; 2]) == [0.09942819 0.29923607; 0.29923607 0.90057181], ρ in :SDP]
-    p = satisfy(constraints)
+    p = satisfy(constraints; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test evaluate(ρ) ≈ [0.09942819 0.29923607 0 0; 0.299237 0.900572 0 0; 0 0 0 0; 0 0 0 0] atol=atol rtol=rtol
@@ -267,11 +299,13 @@ end
     A = randn(m,n) + im*randn(m,n)
     b = A * xo
     x = Variable(n)
-    p1 = minimize(sum(x), A*x == b, x>=0)
+    p1 = minimize(sum(x), A*x == b, x>=0; numeric_type = T)
+
     handle_problem!(p1)
     x1 = x.value
 
-    p2 = minimize(sum(x), real(A)*x == real(b), imag(A)*x==imag(b), x>=0)
+    p2 = minimize(sum(x), real(A)*x == real(b), imag(A)*x==imag(b), x>=0; numeric_type = T)
+
     handle_problem!(p2)
     x2 = x.value
     if test
@@ -286,13 +320,15 @@ end
     A = randn(m,n) + im*randn(m,n)
     b = A * xo
     x = ComplexVariable(n)
-    p1 = minimize(real(sum(x)), A*x == b, real(x)>=0, imag(x)>=0)
+    p1 = minimize(real(sum(x)), A*x == b, real(x)>=0, imag(x)>=0; numeric_type = T)
+
     handle_problem!(p1)
     x1 = x.value
 
     xr = Variable(n)
     xi = Variable(n)
-    p2 = minimize(sum(xr), real(A)*xr-imag(A)*xi == real(b), imag(A)*xr+real(A)*xi == imag(b), xr>=0, xi>=0)
+    p2 = minimize(sum(xr), real(A)*xr-imag(A)*xi == real(b), imag(A)*xr+real(A)*xi == imag(b), xr>=0, xi>=0; numeric_type = T)
+
     handle_problem!(p2)
     #x2 = xr.value + im*xi.value
 
@@ -308,21 +344,23 @@ end
 @add_problem sdp function sdp_Issue_198(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     ρ = HermitianSemidefinite(2)
     constraints = [ρ == [ 1. 0.; 0.  1.]]
-    p = satisfy(constraints)
+    p = satisfy(constraints; numeric_type = T)
+
     handle_problem!(p)
     if test
-        @test p.status == :Optimal
-        @test p.solution.primal ≈ [0.; 1.; 0.; 0.; 1.; zeros(4)] atol=atol rtol=rtol
+        @test p.status == MOI.OPTIMAL
+        @test evaluate(ρ) ≈ [ 1. 0.; 0.  1.] atol=atol rtol=rtol
         @test p.optval ≈ 0 atol=atol rtol=rtol
     end
 end
 
-@add_problem sdp function sdp_norm2_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+@add_problem sdp function sdp_socp_norm2_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     a = 2+4im
     x = ComplexVariable()
     objective = norm2(a-x)
     c1 = real(x)>=0
-    p = minimize(objective,c1)
+    p = minimize(objective,c1; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test p.optval ≈ 0 atol=atol rtol=rtol
@@ -335,12 +373,13 @@ end
     end
 end
 
-@add_problem sdp function sdp_sumsquares_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+@add_problem sdp function sdp_socp_sumsquares_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     a = [2+4im;4+6im]
     x = ComplexVariable(2)
     objective = sumsquares(a-x)
     c1 = real(x)>=0
-    p = minimize(objective,c1)
+    p = minimize(objective,c1; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test p.optval ≈ 0 atol=atol rtol=rtol
@@ -353,12 +392,13 @@ end
     end
 end
 
-@add_problem sdp function sdp_abs_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+@add_problem sdp function sdp_socp_abs_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     a = [5-4im]
     x = ComplexVariable()
     objective = abs(a-x)
     c1 = real(x)>=0
-    p = minimize(objective,c1)
+    p = minimize(objective,c1; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test p.optval ≈ 0 atol=atol rtol=rtol
@@ -378,7 +418,8 @@ end
     x = ComplexVariable(n,n)
     objective = sumsquares(A - x)
     c1 = x in :SDP
-    p = minimize(objective, c1)
+    p = minimize(objective, c1; numeric_type = T)
+
     handle_problem!(p)
     # test that X is approximately equal to posA:
     l,v = eigen(A)
diff --git a/src/problem_depot/problems/sdp_and_exp.jl b/src/problem_depot/problems/sdp_and_exp.jl
index ef11c2a8f..98c74bdea 100644
--- a/src/problem_depot/problems/sdp_and_exp.jl
+++ b/src/problem_depot/problems/sdp_and_exp.jl
@@ -1,6 +1,7 @@
 @add_problem sdp_and_exp function sdp_and_exp_log_det_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2, 2)
-    p = maximize(logdet(x), [x[1, 1] == 1, x[2, 2] == 1])
+    p = maximize(logdet(x), [x[1, 1] == 1, x[2, 2] == 1]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
diff --git a/src/problem_depot/problems/socp.jl b/src/problem_depot/problems/socp.jl
index 5269a6d6e..ef104ae7f 100644
--- a/src/problem_depot/problems/socp.jl
+++ b/src/problem_depot/problems/socp.jl
@@ -1,8 +1,9 @@
-@add_problem socp function socp_norm_2_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+@add_problem socp function socp_dual_norm_2_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2, 1)
     A = [1 2; 2 1; 3 4]
     b = [2; 3; 4]
-    p = minimize(norm2(A * x + b))
+    p = minimize(norm2(A * x + b); numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -16,7 +17,8 @@
     A = [1 2; 2 1; 3 4]
     b = [2; 3; 4]
     lambda = 1
-    p = minimize(norm2(A * x + b) + lambda * norm2(x), x >= 1)
+    p = minimize(norm2(A * x + b) + lambda * norm2(x), x >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -24,11 +26,13 @@
     if test
         @test p.optval ≈ 14.9049 atol=atol rtol=rtol
         @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=atol rtol=rtol
+        @test p.constraints[1].dual ≈ [4.4134, 5.1546] atol=atol rtol=rtol
     end
 
     x = Variable(2)
 
-    p = minimize(norm2([x[1] + 2x[2] + 2; 2x[1] + x[2] + 3; 3x[1]+4x[2] + 4]) + lambda * norm2(x), x >= 1)
+    p = minimize(norm2([x[1] + 2x[2] + 2; 2x[1] + x[2] + 3; 3x[1]+4x[2] + 4]) + lambda * norm2(x), x >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -37,13 +41,15 @@
     if test
         @test p.optval ≈ 14.9049 atol=atol rtol=rtol
         @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=atol rtol=rtol
+        @test p.constraints[1].dual ≈ [4.4134, 5.1546] atol=atol rtol=rtol
     end
 
     x = Variable(2, 1)
     A = [1 2; 2 1; 3 4]
     b = [2; 3; 4]
     lambda = 1
-    p = minimize(norm2(A * x + b) + lambda * norm_1(x), x >= 1)
+    p = minimize(norm2(A * x + b) + lambda * norm_1(x), x >= 1; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -51,13 +57,16 @@
     if test
         @test p.optval ≈ 15.4907 atol=atol rtol=rtol
         @test evaluate(norm2(A * x + b) + lambda * norm_1(x)) ≈ 15.4907 atol=atol rtol=rtol
+        @test p.constraints[1].dual ≈ [4.7062, 5.4475] atol=atol rtol=rtol
+
     end
 end
 
-@add_problem socp function socp_frobenius_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+@add_problem socp function socp_dual_frobenius_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     m = Variable(4, 5)
     c = [m[3, 3] == 4, m >= 1]
-    p = minimize(norm(vec(m), 2), c)
+    p = minimize(norm(vec(m), 2), c; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -65,6 +74,9 @@ end
     if test
         @test p.optval ≈ sqrt(35) atol=atol rtol=rtol
         @test evaluate(norm(vec(m), 2)) ≈ sqrt(35) atol=atol rtol=rtol
+        @test p.constraints[1].dual ≈ 0.6761 atol=atol rtol=rtol
+        dual = 0.1690 .* ones(4, 5); dual[3, 3] = 0
+        @test p.constraints[2].dual ≈ dual atol=atol rtol=rtol
     end
 end
 
@@ -74,7 +86,8 @@ end
     b = [-3; 9; 5]
     c = [3 2 4]
     d = -3
-    p = minimize(quadoverlin(A*x + b, c*x + d))
+    p = minimize(quadoverlin(A*x + b, c*x + d); numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -89,7 +102,8 @@ end
     x = Variable(2, 1)
     A = [1 2; 2 1; 3 4]
     b = [2; 3; 4]
-    p = minimize(sumsquares(A*x + b))
+    p = minimize(sumsquares(A*x + b); numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -104,7 +118,8 @@ end
     x = Variable(2, 1)
     A = [1 2; 2 1; 3 4]
     b = [2; 3; 4]
-    p = minimize(sum(square(A*x + b)))
+    p = minimize(sum(square(A*x + b)); numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -118,7 +133,7 @@ end
     A = [1 2; 2 1; 3 4]
     b = [2; 3; 4]
     expr = A * x + b
-    p = minimize(sum(dot(^)(expr,2))) # elementwise ^
+    p = minimize(sum(dot(^)(expr,2)); numeric_type = T) # elementwise ^
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -128,7 +143,7 @@ end
         @test evaluate(sum(broadcast(^, expr, 2))) ≈ 0.42105 atol=atol rtol=rtol
     end
 
-    p = minimize(sum(dot(*)(expr, expr))) # elementwise *
+    p = minimize(sum(dot(*)(expr, expr)); numeric_type = T) # elementwise *
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -141,7 +156,8 @@ end
 
 @add_problem socp function socp_inv_pos_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(4)
-    p = minimize(sum(invpos(x)), invpos(x) < 2, x > 1, x == 2, 2 == x)
+    p = minimize(sum(invpos(x)), invpos(x) < 2, x > 1, x == 2, 2 == x; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -152,7 +168,8 @@ end
     end
 
     x = Variable(3)
-    p = minimize(sum(dot(/)([3,6,9], x)), x<=3)
+    p = minimize(sum(dot(/)([3,6,9], x)), x<=3; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test x.value ≈ fill(3.0, (3, 1)) atol=atol rtol=rtol
@@ -161,7 +178,8 @@ end
     end
 
     x = Variable()
-    p = minimize(sum([3,6,9]/x), x<=3)
+    p = minimize(sum([3,6,9]/x), x<=3; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test x.value ≈ 3 atol=atol rtol=rtol
@@ -173,18 +191,21 @@ end
 @add_problem socp function socp_geo_mean_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(2)
     y = Variable(2)
-    p = minimize(geomean(x, y), x >= 1, y >= 2)
+    p = minimize(geomean(x, y), x >= 1, y >= 2; numeric_type = T)
+
     # not DCP compliant
     if test
         @test vexity(p) == ConcaveVexity()
     end
-    p = maximize(geomean(x, y), 1 < x, x < 2, y < 2)
+    p = maximize(geomean(x, y), 1 < x, x < 2, y < 2; numeric_type = T)
+
     # Just gave it a vector as an objective, not okay
     if test
         @test_throws Exception handle_problem!(p)
     end
 
-    p = maximize(sum(geomean(x, y)), 1 < x, x < 2, y < 2)
+    p = maximize(sum(geomean(x, y)), 1 < x, x < 2, y < 2; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test p.optval ≈ 4 atol=atol rtol=rtol
@@ -194,13 +215,15 @@ end
 
 @add_problem socp function socp_sqrt_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable()
-    p = maximize(sqrt(x), 1 >= x)
+    p = maximize(sqrt(x), 1 >= x; numeric_type = T)
+
 end
 
 @add_problem socp function socp_quad_form_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3, 1)
     A = [0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422]
-    p = minimize(quadform(x, A), [x >= 1])
+    p = minimize(quadform(x, A), [x >= 1]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -213,7 +236,8 @@ end
     x = Variable(3, 1)
     A = -1.0*[0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422]
     c = [3 2 4]
-    p = maximize(c*x , [quadform(x, A) >= -1])
+    p = maximize(c*x , [quadform(x, A) >= -1]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -226,7 +250,8 @@ end
 
 @add_problem socp function socp_huber_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     x = Variable(3)
-    p = minimize(sum(huber(x, 1)), x >= 2)
+    p = minimize(sum(huber(x, 1)), x >= 2; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -241,7 +266,8 @@ end
     A = [1 2 3; -1 2 3]
     b = A * ones(3)
     x = Variable(3)
-    p = minimize(norm(x, 4.5), [A * x == b])
+    p = minimize(norm(x, 4.5), [A * x == b]; numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -258,7 +284,8 @@ end
     x = Variable(5)
     q = 1.379;  # q norm constraint that generates many inequalities
     qs = q / (q - 1);  # Conjugate to q
-    p = minimize(x' * v)
+    p = minimize(x' * v; numeric_type = T)
+
     p.constraints += (norm(x, q) <= 1)
     if test
         @test vexity(p) == ConvexVexity()
@@ -279,7 +306,8 @@ end
     b = [-1.82041, -1.67516, -0.866884]
     q = 1.5
     xvar = Variable(2)
-    p = minimize(.5 * sumsquares(xvar) + norm(A * xvar - b, q))
+    p = minimize(.5 * sumsquares(xvar) + norm(A * xvar - b, q); numeric_type = T)
+
     if test
         @test vexity(p) == ConvexVexity()
     end
@@ -323,7 +351,8 @@ end
     x = Variable()
     y = Variable()
 
-    p = minimize(x+y, x>=0, y>=0)
+    p = minimize(x+y, x>=0, y>=0; numeric_type = T)
+
     handle_problem!(p)
     if test
         @test p.optval ≈ 0 atol=atol rtol=rtol
@@ -349,7 +378,8 @@ end
     gamma = Variable(Positive())
     fix!(gamma, 0.7)
 
-    p = minimize(norm(x-a) + gamma*norm(x[1:end-1] - x[2:end]))
+    p = minimize(norm(x-a) + gamma*norm(x[1:end-1] - x[2:end]); numeric_type = T)
+
     handle_problem!(p)
     if test
         o1 = p.optval
@@ -376,7 +406,7 @@ end
     x = Variable(2)
     A = [1 2; 2 4];
     b = [3, 6];
-    p = minimize(norm(x, 1), A*x==b)
+    p = minimize(norm(x, 1), A*x==b; numeric_type = T)
 
     if test
         @test vexity(p) == ConvexVexity()
@@ -393,7 +423,7 @@ end
     x = Variable(2)
     A = [1 2; 2 4];
     b = [3, 6];
-    p = minimize(norm(x, 2), A*x==b)
+    p = minimize(norm(x, 2), A*x==b; numeric_type = T)
     
     test && @test vexity(p) == ConvexVexity()
     
@@ -409,7 +439,7 @@ end
     x = Variable(2)
     A = [1 2; 2 4];
     b = [3, 6];
-    p = minimize(norm(x, Inf), A*x==b)
+    p = minimize(norm(x, Inf), A*x==b; numeric_type = T)
 
     test && @test vexity(p) == ConvexVexity()
 
diff --git a/src/problems.jl b/src/problems.jl
index dabddf90e..7fe074342 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -1,81 +1,41 @@
-import MathProgBase
+export Problem, minimize, maximize, satisfy, add_constraint!, add_constraints!
 
-export Problem, Solution, minimize, maximize, satisfy, add_constraint!, add_constraints!
-export Float64OrNothing
-
-const Float64OrNothing = Union{Float64, Nothing}
-
-# TODO: Cleanup
-mutable struct Solution{T<:Number}
-    primal::Array{T, 1}
-    dual::Array{T, 1}
-    status::Symbol
-    optval::T
-    has_dual::Bool
-end
-
-Solution(x::Array{T, 1}, status::Symbol, optval::T) where {T} =
-    Solution(x, T[], status, optval, false)
-Solution(x::Array{T, 1}, y::Array{T, 1}, status::Symbol, optval::T) where {T} =
-    Solution(x, y, status, optval, true)
-
-mutable struct Problem
+mutable struct Problem{T<:Real}
     head::Symbol
     objective::AbstractExpr
     constraints::Array{Constraint}
-    status::Symbol
-    optval::Float64OrNothing
-    model::Union{MathProgBase.AbstractConicModel, Nothing}
-    solution::Solution
+    status::MOI.TerminationStatusCode
+    model::Union{MOI.ModelLike, Nothing}
 
-    function Problem(head::Symbol, objective::AbstractExpr,
-                     model::Union{MathProgBase.AbstractConicModel, Nothing},
-                     constraints::Array=Constraint[])
+    function Problem{T}(head::Symbol, objective::AbstractExpr,
+                     constraints::Array=Constraint[]) where {T <: Real}
         if sign(objective)== Convex.ComplexSign()
-            error("Objective can not be a complex expression")
+            error("Objective cannot be a complex expression")
         else
-            return new(head, objective, constraints, Symbol("not yet solved"), nothing, model)
+            return new(head, objective, constraints, MOI.OPTIMIZE_NOT_CALLED, nothing)
         end
     end
 end
 
-# constructor if model is not specified
-function Problem(head::Symbol, objective::AbstractExpr, constraints::Array=Constraint[],
-                 solver::Union{MathProgBase.AbstractMathProgSolver, Nothing}=nothing)
-    model = solver !== nothing ? MathProgBase.ConicModel(solver) : solver
-    Problem(head, objective, model, constraints)
-end
-
-# If the problem constructed is of the form Ax=b where A is m x n
-# returns:
-# index: n
-# constr_size: m
-# var_to_ranges a dictionary mapping from variable id to (start_index, end_index)
-# where start_index and end_index are the start and end indexes of the variable in A
-function find_variable_ranges(constraints, id_to_variables)
-    index = 0
-    constr_size = 0
-    var_to_ranges = Dict{UInt64, Tuple{Int, Int}}()
-    for constraint in constraints
-        for i = 1:length(constraint.objs)
-            for (id, val) in constraint.objs[i]
-                if !haskey(var_to_ranges, id) && id != objectid(:constant)
-                    var = id_to_variables[id]
-                    if var.sign == ComplexSign()
-                        var_to_ranges[id] = (index + 1, index + 2*length(var))
-                        index += 2*length(var)
-                    else
-                        var_to_ranges[id] = (index + 1, index + length(var))
-                        index += length(var)
-                    end
-                end
-            end
-            constr_size += constraint.sizes[i]
+function Base.getproperty(p::Problem, s::Symbol)
+    if s === :optval
+        if getfield(p, :status) == MOI.OPTIMIZE_NOT_CALLED
+            return nothing
+        else
+            return MOI.get(p.model, MOI.ObjectiveValue())
         end
+    else
+        return getfield(p, s)
     end
-    return index, constr_size, var_to_ranges
 end
 
+dual_status(p::Problem) = MOI.get(p.model, MOI.DualStatus())
+primal_status(p::Problem) = MOI.get(p.model, MOI.PrimalStatus())
+termination_status(p::Problem) = MOI.get(p.model, MOI.TerminationStatus())
+objective_value(p::Problem) = MOI.get(p.model, MOI.ObjectiveValue())
+
+Problem(args...) = Problem{Float64}(args...)
+
 function vexity(p::Problem)
     bad_vex = [ConcaveVexity, NotDcp]
 
@@ -107,121 +67,34 @@ function conic_form!(p::Problem, unique_conic_forms::UniqueConicForms)
     return objective, objective_var.id_hash
 end
 
-function conic_problem(p::Problem)
-    if length(p.objective) != 1
-        error("Objective must be a scalar")
-    end
-
-    # conic problems have the form
-    # minimize c'*x
-    # st       b - Ax \in cones
-    # our job is to take the conic forms of the objective and constraints
-    # and convert them into vectors b and c and a matrix A
-    # one chunk of rows in b and in A corresponds to each constraint,
-    # and one chunk of columns in b and A corresponds to each variable,
-    # with the size of the chunk determined by the size of the constraint or of the variable
-
-    # A map to hold unique constraints. Each constraint is keyed by a symbol
-    # of which atom generated the constraints, and a integer hash of the child
-    # expressions used by the atom
-    unique_conic_forms = UniqueConicForms()
-    objective, objective_var_id = conic_form!(p, unique_conic_forms)
-    constraints = unique_conic_forms.constr_list
-    conic_constr_to_constr = unique_conic_forms.conic_constr_to_constr
-    id_to_variables = unique_conic_forms.id_to_variables
-
-    # var_to_ranges maps from variable id to the (start_index, stop_index) pairs of the columns of A corresponding to that variable
-    # var_size is the sum of the lengths of all variables in the problem
-    # constr_size is the sum of the lengths of all constraints in the problem
-    var_size, constr_size, var_to_ranges = find_variable_ranges(constraints, id_to_variables)
-    c = spzeros(var_size, 1)
-    objective_range = var_to_ranges[objective_var_id]
-    c[objective_range[1]:objective_range[2]] .= 1
-
-    # slot in all of the coefficients in the conic forms into A and b
-    A = spzeros(constr_size, var_size)
-    b = spzeros(constr_size, 1)
-    cones = Tuple{Symbol, UnitRange{Int}}[]
-    constr_index = 0
-    for constraint in constraints
-        total_constraint_size = 0
-        for i = 1:length(constraint.objs)
-            sz = constraint.sizes[i]
-            for (id, val) in constraint.objs[i]
-                if id == objectid(:constant)
-                    for l in 1:sz
-                        b[constr_index + l] = val[1][l] == 0 ? val[2][l] : val[1][l]
-                    end
-                    #b[constr_index + sz + 1 : constr_index + 2*sz] = val[2]
-                else
-                    var_range = var_to_ranges[id]
-                    if id_to_variables[id].sign == ComplexSign()
-                        A[constr_index + 1 : constr_index + sz, var_range[1] : var_range[1] + length(id_to_variables[id])-1] = -val[1]
-                        A[constr_index + 1 : constr_index + sz, var_range[1] + length(id_to_variables[id]) : var_range[2]] = -val[2]
-                    else
-                        A[constr_index + 1 : constr_index + sz, var_range[1] : var_range[2]] = -val[1]
-                    end
-                end
-            end
-            constr_index += sz
-            total_constraint_size += sz
-        end
-        push!(cones, (constraint.cone, constr_index - total_constraint_size + 1 : constr_index))
-    end
-
-    # find integral and boolean variables
-    vartypes = fill(:Cont, length(c))
-    for var_id in keys(var_to_ranges)
-        variable = id_to_variables[var_id]
-        if :Int in variable.sets
-            startidx, endidx = var_to_ranges[var_id]
-            for idx in startidx:endidx
-                vartypes[idx] = :Int
-            end
-        end
-        if :Bin in variable.sets
-            startidx, endidx = var_to_ranges[var_id]
-            for idx in startidx:endidx
-                vartypes[idx] = :Bin
-            end
-        end
-    end
-
-    if p.head == :maximize
-        c = -c
-    end
-
-    return c, A, b, cones, var_to_ranges, vartypes, constraints, id_to_variables, conic_constr_to_constr
-end
-
-Problem(head::Symbol, objective::AbstractExpr, constraints::Constraint...) =
-    Problem(head, objective, [constraints...])
+Problem{T}(head::Symbol, objective::AbstractExpr, constraints::Constraint...) where {T<:Real} =
+    Problem{T}(head, objective, [constraints...])
 
 # Allow users to simply type minimize
-minimize(objective::AbstractExpr, constraints::Constraint...) =
-    Problem(:minimize, objective, collect(constraints))
-minimize(objective::AbstractExpr, constraints::Array{<:Constraint}=Constraint[]) =
-    Problem(:minimize, objective, constraints)
-minimize(objective::Value, constraints::Constraint...) =
-    minimize(convert(AbstractExpr, objective), collect(constraints))
-minimize(objective::Value, constraints::Array{<:Constraint}=Constraint[]) =
-    minimize(convert(AbstractExpr, objective), constraints)
+minimize(objective::AbstractExpr, constraints::Constraint...; numeric_type = Float64) =
+    Problem{numeric_type}(:minimize, objective, collect(constraints))
+minimize(objective::AbstractExpr, constraints::Array{<:Constraint}=Constraint[]; numeric_type = Float64) =
+    Problem{numeric_type}(:minimize, objective, constraints)
+minimize(objective::Value, constraints::Constraint...; numeric_type = Float64) =
+    minimize(convert(AbstractExpr, objective), collect(constraints); numeric_type = numeric_type)
+minimize(objective::Value, constraints::Array{<:Constraint}=Constraint[]; numeric_type = Float64) =
+    minimize(convert(AbstractExpr, objective), constraints; numeric_type = numeric_type)
 
 # Allow users to simply type maximize
-maximize(objective::AbstractExpr, constraints::Constraint...) =
-    Problem(:maximize, objective, collect(constraints))
-maximize(objective::AbstractExpr, constraints::Array{<:Constraint}=Constraint[]) =
-    Problem(:maximize, objective, constraints)
-maximize(objective::Value, constraints::Constraint...) =
-    maximize(convert(AbstractExpr, objective), collect(constraints))
-maximize(objective::Value, constraints::Array{<:Constraint}=Constraint[]) =
-    maximize(convert(AbstractExpr, objective), constraints)
+maximize(objective::AbstractExpr, constraints::Constraint...; numeric_type = Float64) =
+    Problem{numeric_type}(:maximize, objective, collect(constraints))
+maximize(objective::AbstractExpr, constraints::Array{<:Constraint}=Constraint[]; numeric_type = Float64) =
+    Problem{numeric_type}(:maximize, objective, constraints)
+maximize(objective::Value, constraints::Constraint...; numeric_type = Float64) =
+    maximize(convert(AbstractExpr, objective), collect(constraints); numeric_type = numeric_type)
+maximize(objective::Value, constraints::Array{<:Constraint}=Constraint[]; numeric_type = Float64) =
+    maximize(convert(AbstractExpr, objective), constraints; numeric_type = numeric_type)
 
 # Allow users to simply type satisfy (if there is no objective)
-satisfy(constraints::Constraint...) = Problem(:minimize, Constant(0), [constraints...])
-satisfy(constraints::Array{<:Constraint}=Constraint[]) =
-    Problem(:minimize, Constant(0), constraints)
-satisfy(constraint::Constraint) = satisfy([constraint])
+satisfy(constraints::Constraint...; numeric_type = Float64) = Problem{numeric_type}(:minimize, Constant(0), [constraints...])
+satisfy(constraints::Array{<:Constraint}=Constraint[]; numeric_type = Float64) =
+    Problem{numeric_type}(:minimize, Constant(0), constraints)
+satisfy(constraint::Constraint; numeric_type = Float64) = satisfy([constraint]; numeric_type = numeric_type)
 
 # +(constraints, constraints) is defined in constraints.jl
 add_constraints!(p::Problem, constraints::Array{<:Constraint}) =
diff --git a/src/solution.jl b/src/solution.jl
index 30e024688..5e746674b 100644
--- a/src/solution.jl
+++ b/src/solution.jl
@@ -1,214 +1,321 @@
-import MathProgBase
 export solve!
 
-# semantics:
-# * we *load* problem data from Convex into the MathProgBase model
-# * we *populate* objects in Convex with their optimal values
+# Convert from sets used within Convex to MOI sets
+function get_MOI_set(cone, length_inds)
+    if cone == :SDP
+        set = MOI.PositiveSemidefiniteConeTriangle(Int(sqrt(.25 + 2 * length_inds) - .5))
+    elseif cone == :Zero
+        set = MOI.Zeros(length_inds)
+    elseif cone == :Free
+        set = MOI.Reals(length_inds)
+    elseif cone == :NonNeg
+        set = MOI.Nonnegatives(length_inds)
+    elseif cone == :NonPos
+        set = MOI.Nonpositives(length_inds)
+    elseif cone == :SOC
+        set = MOI.SecondOrderCone(length_inds)
+    elseif cone == :SOCRotated
+        set = MOI.RotatedSecondOrderCone(length_inds)
+    elseif cone == :ExpPrimal
+        set = MOI.ExponentialCone()
+    elseif cone == :ExpDual
+        set = MOI.DualExponentialCone()
+    else
+        error("Cone $cone not found somehow; please file an issue.")
+    end
+    return set
+end
 
-function solve!(problem::Problem,
-                s::MathProgBase.AbstractMathProgSolver;
-                kwargs...)
-    # TODO: warn about wiping out old model if warmstart=true?
-    problem.model = MathProgBase.ConicModel(s)
-    return solve!(problem; kwargs...)
+# add_terms! turns a matrix or vector `X` into a collection of `VectorAffineTerm`'s, which it adds to a list (the first argument).
+# This replicates the behavior of the relation `A[row_indices, col_indices] = X`, where `X` could be a vector or matrix.
+# Here, `row_indices` corresponds to the output index of a `VectorAffineTerm`
+# Instead of `col_indices`, here we have their wrapped counterparts, `var_indices` (a vector of `MOI.VariableIndex`'s).
+function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, matrix::SparseMatrixCSC, row_indices, var_indices) where T
+    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
+    I, J, V = findnz(matrix)
+    for  k = eachindex(I, J, V)
+        term = MOI.VectorAffineTerm{T}(CIs[I[k],J[k]][1], MOI.ScalarAffineTerm{T}(V[k], var_indices[CIs[I[k],J[k]][2]]))
+        push!(terms, term)
+    end
 end
 
-function solve!(problem::Problem;
-                warmstart=false,
-                check_vexity=true,
-                verbose=true)
+function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, vector::SparseVector, row_indices, var_indices) where T
+    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
+    I, V = findnz(vector)
+    for k = eachindex(I, V)
+        term =  MOI.VectorAffineTerm{T}(CIs[I[k]][1], MOI.ScalarAffineTerm{T}(V[k], var_indices[CIs[I[k]][2]]))
+        push!(terms, term)
+    end
+end
 
-    if problem.model === nothing
-        throw(ArgumentError(
-            "The provided problem hasn't been initialized with a conic model.
-            You can resolve this by passing in a `AbstractMathProgSolver`. For example,
-            ```
-            using ECOS
-            solve!(problem, ECOSSolver())
-            ```"
-        ))
+function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, matrix::Matrix, row_indices, var_indices) where T
+    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
+    for ci = CartesianIndices(size(matrix))
+        term = MOI.VectorAffineTerm{T}(CIs[ci][1], MOI.ScalarAffineTerm{T}(matrix[ci], var_indices[CIs[ci][2]]))
+        push!(terms, term)
     end
+end
 
-    if check_vexity
-        vex = vexity(problem)
+function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, vector::Vector, row_indices, var_indices) where T
+    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
+    for k = 1:length(vector)
+        term = MOI.VectorAffineTerm{T}(CIs[k][1], MOI.ScalarAffineTerm{T}(vector[k], var_indices[CIs[k][2]]))
+        push!(terms, term)
     end
+end
 
-    c, A, b, cones, var_to_ranges, vartypes, conic_constraints, id_to_variables, conic_constr_to_constr = conic_problem(problem)
 
-    # load MPB conic problem
-    m = problem.model
-    load_problem!(m, c, A, b, cones, vartypes)
-    if warmstart
-        set_warmstart!(m, problem, length(c), var_to_ranges, id_to_variables)
-    end
-    # optimize MPB conic problem
-    MathProgBase.optimize!(m)
+# This function generates a MOI set and VectorAffineFunction corresponding to a constraint that
+# Convex represents with a `ConicConstr` object.
+# The constraint is of the form Ax + b ∈ set
+# where `A` is a matrix, `x` a vector of variables, and `b` a constant vector.
+# MOI represents this type of constraint as a VectorAffineFunction, where `A` is represented
+# by a collection of `VectorAffineTerm`'s.
+function make_MOI_constr(conic_constr::ConicConstr, var_to_indices, id_to_variables, T)
+    total_constraint_size = sum(conic_constr.sizes)
+    constr_index = 0
+    # We create a vector of `VectorAffineTerm`'s to build the constraint.
+    terms =  MOI.VectorAffineTerm{T}[]
+    b = spzeros(T, total_constraint_size)
 
-    # populate the status, the primal (and possibly dual) solution
-    # and the primal (and possibly dual) variables with values
-    populate_solution!(m, problem, var_to_ranges, conic_constraints, id_to_variables, conic_constr_to_constr)
-    if problem.status != :Optimal && verbose
-        @warn "Problem status $(problem.status); solution may be inaccurate."
+    for i = 1:length(conic_constr.objs)
+        sz = conic_constr.sizes[i]
+        for (id, val) in conic_constr.objs[i]
+            if id == objectid(:constant)
+                for l in 1:sz
+                    b[constr_index + l] = ifelse(val[1][l] == 0, val[2][l],  val[1][l])
+                end
+            else
+                var_indices = var_to_indices[id]
+                if id_to_variables[id].sign == ComplexSign()
+                    l = length(var_indices) ÷ 2
+                    # We create MOI terms and add them to `terms`.
+                    # Real part:
+                    add_terms!(terms, val[1], constr_index + 1 : constr_index + sz, var_indices[1:l])
+
+                    # Imaginary part:
+                    add_terms!(terms, val[2], constr_index + 1 : constr_index + sz, var_indices[l+1 : end])
+
+                else
+                    # Only a real part:
+                    add_terms!(terms, val[1], constr_index + 1 : constr_index + sz, var_indices)
+                end
+            end
+        end
+        constr_index += sz
     end
+    set = get_MOI_set(conic_constr.cone, total_constraint_size)
+    constr_fn =  MOI.VectorAffineFunction{T}(terms, b)
+    return set, constr_fn
 end
 
-function set_warmstart!(m::MathProgBase.AbstractConicModel,
-                        problem::Problem,
-                        n::Int, # length of primal (conic) solution
-                        var_to_ranges,
-                        id_to_variables)
-     # use previously cached solution, if any,
-     try
-         primal = problem.solution.primal
-     catch
-         @warn "Unable to use cached solution to warmstart problem.
-                (Perhaps this is the first time you're solving this problem?)
-                Warmstart may be ineffective."
-         primal = zeros(n)
-     end
-     if length(primal) != n
-         @warn "Unable to use cached solution to warmstart problem.
-                (Perhaps the number of variables or constraints in the problem have changed since you last solved it?)
-                Warmstart may be ineffective."
-         primal = zeros(n)
-     end
-
-     # grab any variables whose values the user may be trying to set
-     load_primal_solution!(primal, var_to_ranges, id_to_variables)
-
-     # notify the model that we're trying to warmstart
-     try
-         MathProgBase.setwarmstart!(m, primal)
-     catch
-         @warn "Unable to warmstart solution.
-               (Perhaps the solver doesn't support warm starts?)
-               Using a cold start instead."
-     end
-     m
-end
 
-function load_problem!(m::MathProgBase.AbstractConicModel, c, A, b, cones, vartypes)
-    # no conic constraints on variables
-    var_cones = fill((:Free, 1:size(A, 2)),1)
-    MathProgBase.loadproblem!(m, vec(Array(c)), A, vec(Array(b)), cones, var_cones)
-
-    # add integer and binary constraints on variables
-    if !all(==(:Cont), vartypes)
-        try
-            MathProgBase.setvartype!(m, vartypes)
-        catch
-            error("model $(typeof(m)) does not support variables of some of the following types: $(unique(vartypes))")
+# This function gets MOI variable indices for each variable used as a key
+# inside a `ConicObj` used in the problem.
+function get_variable_indices!(model, conic_constraints, id_to_variables)
+    var_to_indices = Dict{UInt64,Vector{MOI.VariableIndex}}()
+    for conic_constr in conic_constraints
+        for i = 1:length(conic_constr.objs)
+            for (id, val) in conic_constr.objs[i]
+                if !haskey(var_to_indices, id) && id != objectid(:constant)
+                    var = id_to_variables[id]
+                    if var.sign == ComplexSign()
+                        var_size = 2 * length(var)
+                    else
+                        var_size = length(var)
+                    end
+                    var_to_indices[id] = MOI.add_variables(model, var_size)
+                end
+            end
         end
     end
-    m
+    return var_to_indices
 end
 
-function populate_solution!(m::MathProgBase.AbstractConicModel,
-                            problem::Problem,
-                            var_to_ranges,
-                            conic_constraints,
-                            id_to_variables,
-                            conic_constr_to_constr)
-    dual = try
-        MathProgBase.getdual(m)
-    catch
-        fill(NaN, MathProgBase.numconstr(m))
+
+
+# This function traverses the problem tree and loads
+# an MOI model with the corresponding problem instance.
+# The problem is written in a simple epigraph form.
+function load_MOI_model!(model, problem::Problem{T}) where {T}
+    if length(problem.objective) != 1
+        error("Objective must be a scalar")
     end
 
-    solution = try
-        MathProgBase.getsolution(m)
-    catch
-        fill(NaN, MathProgBase.numvar(m))
+    unique_conic_forms = UniqueConicForms()
+    objective, objective_var_id = conic_form!(problem, unique_conic_forms)
+    conic_constraints = unique_conic_forms.constr_list
+    conic_constr_to_constr = unique_conic_forms.conic_constr_to_constr
+    id_to_variables = unique_conic_forms.id_to_variables
+
+    # `var_to_indices` maps from variable id to the MOI `VariableIndex`'s corresponding to the variable
+    var_to_indices = get_variable_indices!(model, conic_constraints, id_to_variables)
+
+    # the objective: maximize or minimize a scalar variable
+    objective_index = var_to_indices[objective_var_id][] # get the `MOI.VariableIndex` corresponding to the objective
+    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(objective_index))
+    MOI.set(model, MOI.ObjectiveSense(), problem.head == :maximize ? MOI.MAX_SENSE : MOI.MIN_SENSE)
+
+    # Constraints: Generate a MOI function and a MOI sets for each `ConicConstr` object in the problem
+    MOI_constr_fn = Union{MOI.VectorAffineFunction{T},MOI.SingleVariable}[]
+    MOI_sets = Any[]
+    for conic_constr in conic_constraints
+        set, constr_fn = make_MOI_constr(conic_constr, var_to_indices, id_to_variables, T)
+        push!(MOI_sets, set)
+        push!(MOI_constr_fn, constr_fn)
     end
 
-    objective = try
-        MathProgBase.getobjval(m)
-    catch
-        NaN
+    # Add integral and boolean constraints
+    for var_id in keys(var_to_indices)
+        variable = id_to_variables[var_id]
+        if :Int in variable.sets
+            var_indices = var_to_indices[var_id]
+            for idx = eachindex(var_indices)
+                push!(MOI_constr_fn, MOI.SingleVariable(var_indices[idx]))
+                push!(MOI_sets, MOI.Integer())
+            end
+        end
+        if :Bin in variable.sets
+            var_indices = var_to_indices[var_id]
+            for idx in eachindex(var_indices)
+                push!(MOI_constr_fn, MOI.SingleVariable(var_indices[idx]))
+                push!(MOI_sets, MOI.ZeroOne())
+            end
+        end
     end
 
-    if any(isnan, dual)
-        problem.solution = Solution(solution, MathProgBase.status(m), objective)
-    else
-        problem.solution = Solution(solution, dual, MathProgBase.status(m), objective)
+    # Add all the constraints to the model and collect the corresponding MOI indices
+    constraint_indices = MOI.add_constraints(model, MOI_constr_fn, MOI_sets)
+
+    return id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices
+end
+
+
+function solve!(problem::Problem{T}, optimizer::MOI.ModelLike;
+    check_vexity = true,
+    verbose = true,
+    warmstart = false) where {T}
+
+    if check_vexity
+        vex = vexity(problem)
     end
 
-    populate_variables!(problem, var_to_ranges, id_to_variables)
+    model = MOIB.full_bridge_optimizer(
+        MOIU.CachingOptimizer(
+            MOIU.UniversalFallback(MOIU.Model{T}()),
+            optimizer
+        ),
+        T
+    )
+    
+    id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices = load_MOI_model!(model, problem)
 
-    if problem.solution.has_dual
-        populate_duals!(conic_constraints, problem.solution.dual, conic_constr_to_constr)
+    if warmstart
+        warmstart_variables!(model, var_to_indices, id_to_variables, T, verbose)
     end
+ 
+    MOI.optimize!(model)
+    problem.model = model
 
-    # minimize -> maximize
-    if (problem.head == :maximize)
-        problem.solution.optval = -problem.solution.optval
+    # populate the status, primal variables, and dual variables (when possible)
+    moi_populate_solution!(model, problem, id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices)
+    
+    if problem.status != MOI.OPTIMAL && verbose
+        @warn "Problem status $(problem.status); solution may be inaccurate."
     end
+end
 
-    # Populate the problem with the solution
-    problem.optval = problem.solution.optval
-    problem.status = problem.solution.status
 
-    problem
+function warmstart_variables!(model, var_to_indices, id_to_variables, T, verbose)
+    if MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
+        for (id, var_inds) in pairs(var_to_indices)
+            x = id_to_variables[id]
+            value = x.value
+            value === nothing && continue
+            value_vec = packvec(value, sign(x) == ComplexSign())
+            value_vec = convert(Vector{T}, value_vec)
+            MOI.set(model, MOI.VariablePrimalStart(), var_inds, value_vec)
+        end
+    elseif verbose
+        @warn "Skipping variable warmstart; the solver does not support it."
+    end
 end
 
-function populate_variables!(problem::Problem, var_to_ranges::Dict{UInt64, Tuple{Int, Int}}, id_to_variables)
-    x = problem.solution.primal
-    for (id, (start_index, end_index)) in var_to_ranges
-        var = id_to_variables[id]
-        sz = var.size
-        if var.sign != ComplexSign()
-            var.value = reshape(x[start_index:end_index], sz[1], sz[2])
-            if sz == (1, 1)
-                var.value = var.value[1]
-            end
-        else
-            real_value = reshape(x[start_index:start_index + div(end_index-start_index+1,2)-1], sz[1], sz[2])
-            imag_value = reshape(x[start_index + div(end_index-start_index+1,2):end_index], sz[1], sz[2])
-            var.value = real_value + im*imag_value
-            if sz == (1, 1)
-                var.value = var.value[1]
-            end
+function moi_populate_solution!(model::MOI.ModelLike, problem, id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices)
+    status = MOI.get(model, MOI.TerminationStatus())
+    problem.status = status
+
+    ## This can throw if the problem was not solved correctly
+    # objective = MOI.get(model, MOI.ObjectiveValue())
+    # problem.optval = objective
+    ## Instead, we use a `getproperty` override to access the objective value only when the user asks.
+    ## That way `moi_populate_solution!` will complete successfully.
+
+    dual_status = MOI.get(model, MOI.DualStatus())
+    primal_status = MOI.get(model, MOI.PrimalStatus())
+
+    if primal_status != MOI.NO_SOLUTION
+        for (id, var_indices) in var_to_indices
+            var = id_to_variables[id]
+            vectorized_value =  MOI.get(model, MOI.VariablePrimal(), var_indices)
+            var.value = unpackvec(vectorized_value, size(var), var.sign == ComplexSign())
+        end
+    else
+        for (id, var_indices) in var_to_indices
+            var = id_to_variables[id]
+            var.value = nothing
         end
     end
-end
 
-# populates the solution vector from the .value fields of variables
-# for use in warmstarting
-# TODO: it would be super cool to grab the other expressions that appear in the primal solution vector,
-# get their `expression_to_range`,
-# and populate them too using `evaluate`
-function load_primal_solution!(primal::Array{Float64,1}, var_to_ranges::Dict{UInt64, Tuple{Int, Int}}, id_to_variables)
-    for (id, (start_index, end_index)) in var_to_ranges
-        var = id_to_variables[id]
-        if var.value !== nothing
-            sz = size(var.value)
-            if length(sz) <= 1
-                primal[start_index:end_index] = var.value
-            else
-                primal[start_index:end_index] = reshape(var.value, sz[1]*sz[2], 1)
-            end
+    if dual_status != MOI.NO_SOLUTION
+        for (idx, conic_constr) in enumerate(conic_constraints)
+            haskey(conic_constr_to_constr, conic_constr) || continue
+            constr = conic_constr_to_constr[conic_constr]
+            MOI_constr_indices = constraint_indices[idx]
+            dual_value_vectorized = MOI.get(model, MOI.ConstraintDual(), MOI_constr_indices)
+            iscomplex = sign(constr.lhs) == ComplexSign() || sign(constr.rhs) == ComplexSign()
+            constr.dual = unpackvec(dual_value_vectorized, constr.size, iscomplex)
+        end
+    else
+        for (idx, conic_constr) in enumerate(conic_constraints)
+            haskey(conic_constr_to_constr, conic_constr) || continue
+            constr = conic_constr_to_constr[conic_constr]
+            constr.dual = nothing
         end
     end
+
 end
 
-function populate_duals!(constraints::Array{ConicConstr}, dual::Vector, conic_constr_to_constr)
-    constr_index = 1
-    for constraint in constraints
-        # conic_constr_to_constr only has keys for conic constraints with a single objective
-        # so this will work
-        if haskey(conic_constr_to_constr, constraint)
-            sz = constraint.sizes[1]
-            c = conic_constr_to_constr[constraint]
-            c.dual = reshape(dual[constr_index:constr_index+sz-1], c.size)
-            if c.size == (1, 1)
-                c.dual = c.dual[1]
-            end
-            constr_index += sz
-        else
-            for i = 1:length(constraint.objs)
-                sz = constraint.sizes[i]
-                constr_index += sz
-            end
-        end
+
+# This is type unstable!
+function unpackvec(v::AbstractVector, size::Tuple{Int,Int}, iscomplex::Bool)
+    if iscomplex && length(v) == 2
+        return v[1] + im * v[2]
+    elseif iscomplex
+        l = length(v) ÷ 2
+        # use views?
+        return reshape(v[1:l], size) + im * reshape(v[l+1 : end], size)
+    elseif length(v) == 1
+        return v[]
+    else
+        return reshape(v, size)
+    end
+end
+
+# In `packvec`, we use `real` even in the `iscomplex == false` branch
+# for type stability. Note here `iscomplex` refers to the `Sign` of the variable
+# associated to the values here, not whether or not the variable actually holds
+# a complex number at this point. Hence, one could have `iscomplex==true` and
+# `isreal(value)==true` or even `value isa Real` could hold.
+function packvec(value::Number, iscomplex::Bool)
+    iscomplex ? [real(value), imag(value)] : [real(value)]
+end
+
+function packvec(value::AbstractArray, iscomplex::Bool)
+    value = reshape(value, length(value))
+    if iscomplex
+        return [real(value); imag(value)]
+    else
+        return real(value) 
     end
 end
diff --git a/src/solver_info.jl b/src/solver_info.jl
deleted file mode 100644
index f367fe5e7..000000000
--- a/src/solver_info.jl
+++ /dev/null
@@ -1,50 +0,0 @@
-import MathProgBase
-export can_solve_mip, can_solve_socp, can_solve_sdp, can_solve_exp
-
-# NOTE: Convex has been tested with the following solvers:
-#   ECOS:   ECOSSolver
-#   SCS:    SCSSolver
-#   Gurobi: GurobiSolver
-#   Mosek:  MosekSolver
-#   GLPKMathProgInterface: GLPKSolverMIP
-# It has not been tested with other solvers such a CPLEX.
-
-function can_solve_mip(solver)
-    name = typeof(solver).name.name
-    if name == :GurobiSolver || name == :MosekSolver || name == :GLPKSolverMIP || name == :CPLEXSolver || name == :CbcSolver
-        return true
-    else
-        @info "$name cannot solve mixed integer programs. Consider using Gurobi, Mosek, or GLPK."
-        return false
-    end
-end
-
-function can_solve_socp(solver)
-    if :SOC in MathProgBase.supportedcones(solver)
-        return true
-    else
-        name = typeof(solver).name.name
-        @info "$name cannot solve second order cone programs. Consider using SCS, ECOS, Mosek, or Gurobi."
-        return false
-    end
-end
-
-function can_solve_exp(solver)
-    if :ExpPrimal in MathProgBase.supportedcones(solver)
-        return true
-    else
-        name = typeof(solver).name.name
-        @info "$name cannot solve exponential programs. Consider using SCS or ECOS."
-        return false
-    end
-end
-
-function can_solve_sdp(solver)
-    if :SDP in MathProgBase.supportedcones(solver)
-        return true
-    else
-        name = typeof(solver).name.name
-        @info "$name cannot solve semidefinite programs. Consider using SCS or Mosek."
-        return false
-    end
-end
diff --git a/src/utilities/show.jl b/src/utilities/show.jl
index ff0235dfb..fcaa69258 100644
--- a/src/utilities/show.jl
+++ b/src/utilities/show.jl
@@ -184,7 +184,13 @@ end
 
 function show(io::IO, p::Problem)
     TreePrint.print_tree(io, p, MAXDEPTH[], MAXWIDTH[])
-    print(io, "\ncurrent status: $(p.status)")
+    if p.status == MOI.OPTIMIZE_NOT_CALLED
+        print(io, "\nstatus: `solve!` not called yet")
+    else
+        print(io, "\ntermination status: $(p.status)")
+        print(io, "\nprimal status: $(primal_status(p))")
+        print(io, "\ndual status: $(dual_status(p))")
+    end
     if p.status == "solved"
         print(io, " with optimal value of $(round(p.optval, digits=4))")
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0d9e542ad..c9b1dbcff 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,20 +1,26 @@
 using Convex
 using Convex.ProblemDepot: run_tests
 using Test
-using SCS, ECOS, GLPKMathProgInterface
+using SCS, ECOS, GLPK
 
+using MathOptInterface
+const MOI = MathOptInterface
+const MOIU = MOI.Utilities
 
 # Seed random number stream to improve test reliability
 using Random
 Random.seed!(2)
 
 @testset "ProblemDepot" begin
-    @testset "Problems can run without `solve!`ing if `test==false`" begin
+    @testset "Problems can run without `solve!`ing if `test==false`; T=$T" for T in (Float64, BigFloat)
         Convex.ProblemDepot.foreach_problem() do name, func
             @testset "$name" begin
                 # We want to check to make sure this does not throw
-                func(Convex.conic_problem, Val(false), 0.0, 0.0, Float64)
-                @test true
+                func(Val(false), 0.0, 0.0, T) do problem
+                    @test problem isa Convex.Problem{T} # check numeric type
+                    model = MOIU.MockOptimizer(MOIU.Model{T}())
+                    Convex.load_MOI_model!(model, problem) # make sure it loads without throwing
+                end
             end
         end
     end
@@ -23,21 +29,29 @@ end
 @testset "Convex" begin
     include("test_utilities.jl")
 
+    @testset "SCS with warmstarts" begin
+        # We exclude `sdp_matrix_frac_atom` due to the bug https://github.com/JuliaOpt/SCS.jl/issues/153
+        run_tests(; exclude=[r"mip", r"sdp_matrix_frac_atom"]) do p
+            solve!(p, SCS.Optimizer(verbose=0, eps=1e-6); warmstart = true)
+        end
+    end
+
     @testset "SCS" begin
-        run_tests(; exclude=[r"mip"]) do p
-            solve!(p, SCSSolver(verbose=0, eps=1e-6))
+        # We exclude `sdp_matrix_frac_atom` due to the bug https://github.com/JuliaOpt/SCS.jl/issues/153
+        run_tests(; exclude=[r"mip", r"sdp_matrix_frac_atom"]) do p
+            solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
         end
     end
 
     @testset "ECOS" begin
         run_tests(; exclude=[r"mip", r"sdp"]) do p
-            solve!(p, ECOSSolver(verbose=0))
+            solve!(p, ECOS.Optimizer(verbose=0))
         end
     end
 
-    @testset "GLPK MIP" begin
-        run_tests(; exclude=[r"socp", r"sdp", r"exp", r"dual"]) do p
-            solve!(p, GLPKSolverMIP())
+    @testset "GLPK" begin
+        run_tests(; exclude=[r"exp", r"sdp", r"socp"]) do p
+            solve!(p, GLPK.Optimizer(msg_lev = GLPK.OFF))
         end
     end
 end
diff --git a/test/test_utilities.jl b/test/test_utilities.jl
index 28f7f7115..e62f86597 100644
--- a/test/test_utilities.jl
+++ b/test/test_utilities.jl
@@ -1,5 +1,12 @@
 @testset "Utilities" begin
 
+    @testset "`solve!` does not return anything" begin
+        x = Variable()
+        p = satisfy(x >= 0)
+        output = solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
+        @test output === nothing
+    end
+
     @testset "Show" begin
         x = Variable()
         @test sprint(show, x) == """
@@ -37,7 +44,7 @@
            ├─ real variable ($(Convex.show_id(x)))
            └─ 3
         
-        current status: not yet solved"""
+        status: `solve!` not called yet"""
         
         x = ComplexVariable(2,3)
         @test sprint(show, x) == """
@@ -48,22 +55,79 @@
         $(Convex.show_id(x))"""
 
         # test `MAXDEPTH`
+        # We construct a binary tree of depth >= 3
+        # to make sure it gets truncated appropriately.
         x = Variable(2)
         y = Variable(2)
-        p = minimize(sum(x), hcat(hcat(hcat(hcat(x,y), hcat(x,y)),hcat(hcat(x,y), hcat(x,y))),hcat(hcat(hcat(x,y), hcat(x,y)),hcat(hcat(x,y), hcat(x,y)))) == hcat(hcat(hcat(hcat(x,y), hcat(x,y)),hcat(hcat(x,y), hcat(x,y))),hcat(hcat(hcat(x,y), hcat(x,y)),hcat(hcat(x,y), hcat(x,y)))))
-        @test sprint(show, p) == "minimize\n└─ sum (affine; real)\n   └─ 2-element real variable ($(Convex.show_id(x)))\nsubject to\n└─ == constraint (affine)\n   ├─ hcat (affine; real)\n   │  ├─ hcat (affine; real)\n   │  │  ├─ …\n   │  │  └─ …\n   │  └─ hcat (affine; real)\n   │     ├─ …\n   │     └─ …\n   └─ hcat (affine; real)\n      ├─ hcat (affine; real)\n      │  ├─ …\n      │  └─ …\n      └─ hcat (affine; real)\n         ├─ …\n         └─ …\n\ncurrent status: not yet solved" 
+        level3 = hcat(x,y)
+        level2 = hcat(level3, level3)
+        root = hcat(level2, level2)
+        p = minimize(sum(x), root == root)
+        @test sprint(show, p) == """
+        minimize
+        └─ sum (affine; real)
+           └─ 2-element real variable ($(Convex.show_id(x)))
+        subject to
+        └─ == constraint (affine)
+           ├─ hcat (affine; real)
+           │  ├─ hcat (affine; real)
+           │  │  ├─ …
+           │  │  └─ …
+           │  └─ hcat (affine; real)
+           │     ├─ …
+           │     └─ …
+           └─ hcat (affine; real)
+              ├─ hcat (affine; real)
+              │  ├─ …
+              │  └─ …
+              └─ hcat (affine; real)
+                 ├─ …
+                 └─ …
+
+        status: `solve!` not called yet"""
 
         # test `MAXWIDTH`
         x = Variable()
         p = satisfy([ x == i for i = 1:100])
-        @test sprint(show, p) == "minimize\n└─ 0\nsubject to\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 1\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 2\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 3\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 4\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 5\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 6\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 7\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 8\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 9\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 10\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 11\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 12\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 13\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 14\n├─ == constraint (affine)\n│  ├─ real variable ($(Convex.show_id(x)))\n│  └─ 15\n⋮\n\ncurrent status: not yet solved"
+        old_maxwidth = Convex.MAXWIDTH[]
+        Convex.MAXWIDTH[] = 2
+        @test sprint(show, p) == """
+            minimize
+            └─ 0
+            subject to
+            ├─ == constraint (affine)
+            │  ├─ real variable ($(Convex.show_id(x)))
+            │  └─ 1
+            ├─ == constraint (affine)
+            │  ├─ real variable ($(Convex.show_id(x)))
+            │  └─ 2
+            ⋮
+
+            status: `solve!` not called yet"""
+        Convex.MAXWIDTH[] = old_maxwidth
+
+        # solved problem
+        x = Variable()
+        p = satisfy(x >= 0)
+        output = solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
+        @test sprint(show, p) == """
+                minimize
+                └─ 0
+                subject to
+                └─ >= constraint (affine)
+                   ├─ real variable ($(Convex.show_id(x)))
+                   └─ 0
+
+                termination status: OPTIMAL
+                primal status: FEASIBLE_POINT
+                dual status: FEASIBLE_POINT"""
     end
 
     @testset "clearmemory" begin
         @test_deprecated Convex.clearmemory()
     end
 
-    @testset "ConicObj" for T = [UInt32, UInt64]
+    @testset "ConicObj with type $T" for T = [UInt32, UInt64]
         c = ConicObj()
         z = zero(T)
         @test !haskey(c, z)